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LSPRAY-IV: A Lagrangian Spray Module 


M.S. Raju 
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Cleveland, Ohio 44135 


ABSTRACT 

LSPRAY-IV is a Lagrangian spray solver devel- 
oped for application with parallel computing and un- 
structured grids. It is designed to be massively paral- 
lel and could easily be coupled with any existing gas- 
phase flow and/or Monte Carlo Probability Density 
Function (PDF) solvers. The solver accommodates 
the use of an unstructured mesh with mixed elements 
of either triangular, quadrilateral, and/or tetrahedral 
type for the gas flow grid representation. It is mainly 
designed to predict the flow, thermal and transport 
properties of a rapidly vaporizing spray because of 
its importance in aerospace application. The manual 
provides the user with an understanding of various 
models involved in the spray formulation, its code 
structure and solution algorithm, and various other 
issues related to parallelization and its coupling with 
other solvers. With the development of LSPRAY-IV, 
we have advanced the state-of-the-art in spray com- 
putations in several important ways. Some of the 
major features of the spray module are: 

• It facilitates the use of both unstructured grids 
and parallel computing and, thereby, facilitates 
large-scale combustor computations involving 
complex geometrical configurations. 

• In order to deal with modern gas-turbine fuels 
that are mixtures of many compounds, we have 
extended the spray formulation to the modeling 
of multicomponent liquid fuels. 

• In order to extend the applicability of the spray 
computations over a wide range of low-pressure 
conditions, we have completed the implemen- 
tation of the liquid-phase variable thermo- 
transport properties into our spray formulation. 

• In order to reduce some uncertainty associated 
with the specification of the initial droplet con- 
ditions, we have completed the integration of an 


atomization module into our spray solution pro- 
cedure. The atomization module contains the 
following primary atomization models: (a) sheet 
breakup, (b) air blast, (c) blob jet, and (d) BLS 
(Boundary-Layer Stripping), and the following 
secondary droplet breakup models: (a) Rayleigh- 
Taylor, (b) TAB (Taylor Analogy Breakup), and 
(c) ETAB (Enhanced Taylor Analogy Breakup). 

• Because of its importance in some NASA-related 
critical applications, we have extended the solu- 
tion procedure to the modeling of superheated 
sprays. 

• The spray module is designed in such a way so 
that it could easily be coupled with other CFD 
codes. 

• It can be used in the calculation of both steady 
as well as unsteady computations. 

• We have developed and implemented a time- 
averaging method into the calculation of the 
liquid-phase source contributions in order to ac- 
celerate convergence towards a steady state so- 
lution. 

• The spray module has a multi-liquid and multi- 
injection capability. 

With the aim of improving the overall solution 
procedure involving sprays, we have made several rel- 
evant contributions to the gas side of the computa- 
tions: 

• In order to demonstrate the importance of chem- 
istry/turbulence interactions in the modeling 
of reacting sprays, we have extended the joint 
scalar Monte Carlo PDF (Probability Density 
Function) approach to the modeling of spray 
flames, unstructured grids, and parallel comput- 
ing. 
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• In order to account for the nonideal gas behav- 
ior under critical and supercritical conditions, 
we have completed the integration of a high- 
pressure EOS (Equation-of-State) into our CFD 
module. Also, we have completed the implemen- 
tation of a high-pressure correction into the cal- 
culation of the gas-phase transport properties. 

The spray solution procedure provided favor- 
able results when applied to the modeling of sev- 
eral different spray/gaseous flames - representa- 
tive of those encountered in gas-turbine combus- 
tors, stratified-charge rotary combustion (Wankel) 
engines, supersonic diffusion flames, and pulse det- 
onation combustion devices. Its use has been demon- 
strated in various important NASA projects: the 
NASA’s fundamental aeronautics program initiative 
on emissions through the Supersonics (SUP) and 
Subsonic Fixed Wing (SFW) project offices, Ultra- 
Efficient Engine Technology (UEET), Pulse Deto- 
nation Combustion Technology (PDCT), & Rotary 
Combustion Engine Technology Enablement Project 
(RCETEP). 
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NOMENCLATURE 


a 

liquid jet radius, m, or 

m 

rh k 

rilk, flash 

CL 0 

parent drop radius, m 
initial droplet radius, m 

rilko 


outward area normal vector 


of the nth surface, m 2 

rh k ,t 

A 

a constant 

Bk 

Spalding mass transfer number 

m Q 

B t 

Spalding heat transfer number 

m(t) 

B 0 

a constant 

M 

Cfreq 

a constant 

M a 

c 

a constant 


c D 

drag coefficient 

M f 

Mi 

c p 

specific heat, J / (kg K) 

D 

turbulent diffusion coefficient, m 2 /s 

Dab 

diffusion coefficient, cm 2 /s 

M sbe d 

d 

droplet diameter, m 

MMD 

db 

droplet diameter after primary breakup, m 

Tlk 

dL 

ligament diameter, m 

"'parent 

do 

orifice diameter, m 

product 

n(t) 

dparent. 

parent drop diameter, m 

dproduct 

product drop diameter, m 

N 

dt 

time increment, s 

N n 

d6 

half-cone angle, deg. 

Nt 

Dparent 

energy contained in the parent 

J 


drop, (kg m 2 )/s 2 

N v 

Dproduct 

energy contained in the product 

P 

Nu 


drop, (kg m 2 )/s 2 

P 

h 

specific enthalpy, J/kg, 

Pc 


liquid sheet thickness, m 

Pr 

loi h 

modified Bessel functions of the first kind 

Pr 

k 

turbulence kinetic energy, m 2 /s 2 , 

Psat 


wave number, 1/m, or 

Q 


thermal conductivity, J/(ms K) 

V 

kij 

binary interaction parameter 

Tk 

Vko 

K b 

a breakup constant, 1/s 

K sb 

wave number associated with 

sheet breakup, 1/m 

tsmr 

/ \ 

ki,k 2 

constants 

^ U 

Re 

K Lb 

wave number associated with 

RND 


ligament breakup, 1/m 

Sh 

Dq , K { 

modified Bessel functions of the 

second kind 

Sk 

l 

modified wave number, 1/m 

Smlc 

Ik 

mixture latent heat of evaporation, J /kg 

Smle 

lk,eff 

effective latent heat of evaporation, 


Ikn 

J/kg (defined in Eq. (35)) 
heat of vaporization at normal 

Smlm 

L 

boiling point, J/kg 

liquid sheet breakup length, m 



liquid mass flow rate, kg/s 

droplet vaporization rate, kg/s 

droplet vaporization rate under 

flash evaporating conditions, kg/s 

initial mass flow rate associated 

with kth droplet group, kg/s 

droplet vaporization rate due 

to heat transfer, kg/s 

mass of the parent drop, kg 

mean mass of the product drop distribution, kg 

molecular weight, kg/kg-mole 

molecular weight of gas excluding 

fuel species, kg/kg-mole 

molecular weight of fuel, kg/kg-mole 

molecular weight of ith 

species, kg/kg-mole 

shed drop mass, kg 

mass mean diameter, m 

number of droplets in kth group 

drop number in the parent group 

drop number in the product group 

non-dimensional number (=mo/fh(i)) 

drop number 

parent drop number 

number of surfaces contained in 

a given computational cell 

total number of computational cells 

Nusselt number 

pressure, N/m 2 

critical pressure, N/m 2 

Prandtl number 

reduced pressure, P/P c 

saturation pressure, N/m 2 

density ratio {—p g /pi) 

radius (=- v /a 3 /tsmr), m 

droplet radius, m 

initial droplet radial location, m 

Sauter mean radius, m 

gas constant, J/(kg K) 

Reynolds number 

random number 

Sherwood number 

droplet radius-squared ( = r 2 .), m 2 

source term contribution from liquid 

exchange to mass conservation, kg/s 

source term contribution from liquid 

exchange to energy conservation, J/s 

source term contribution from liquid 

exchange to momentum conservation, kg m/s 2 
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Smis source term contribution from liquid 
exchange to species conservation, kg/s 
SMD Sauter mean diameter, m 
t time, s 

td non-dimensional time (= 2p /° ) 

T temperature, K, or 

non-dimensional number (=ZWe g ) 

Tf, boiling temperature, K 

T n b normal boiling temperature, K 

T c critical temperature, K 

7i. kth droplet temperature, K 

T r reduced temperature, T /T c 

iid drop deceleration, m/s 2 

Ui ith velocity component, m/s 

Uik ith velocity component of 

kth droplet group, m/s 

U relative velocity between gas and liquid, 
m/s, or relative drop velocity, m/s 
vt product drop velocity component, m/s 

V initial liquid velocity, m/s 

V c volume of the computational cell, m 3 , or 

critical molar volume, m 3 /kg-mole 
V n molar volume at normal pressure, 

m 3 /kg- mole 

We Weber number (— Pa ^ f/ ) 

Wj gas-phase chemical reaction rate, 1 /s 
x Cartesian coordinate, m, or drop 

deformation distance, m 

Xi Cartesian coordinate in the ith direction, m 

y Cartesian coordinate, m, or non-dimensional 

deformation parameter (= — ) 
y : i mass fraction of jth species 

x spatial vector 

z Cartesian coordinate, m 

Z compressibility factor (= 7 ^), or 

non-dimensional number (= ) 

Z c critical compressibility factor (=Sdk) 

a represents a coordinate related to 

a Hill’s Vortex, spray cone rotation 
angle, deg., or a constant 

a s overall heat transfer coefficient, kJ/sm 2o K 

a constant 
«2 a constant 

/3 spray cone rotation angle, deg. 

X mole fraction 

5 Dirac-delta function, or 

initial liquid sheet thickness, m 
A p pressure drop in the injector, N/m 2 

A tf local time step used in the flow solver, s 


A t g i global time step used in the spray solver, s 

A tu fuel injection time step, s 

At m i allowable time step in the spray solver, s 
AV computational cell volume, m 3 

e rate of turbulence dissipation, m 2 /s 3 

ej fractional mass evaporating rate of species 
at the droplet surface 
77 wave amplitude, m 

turbulent diffusion coefficient, kg/ms, or 

a factor (=210 [ Tc pl ] ^ 6 ) 

A thermal conductivity, J / (ms K) , or 
wavelength, m 

A wavelength associated with fi, or 

wavelength associated with 
Rayleigh-Taylor breakup, m 
/i dynamic viscosity, kg/ms 

v kinematic viscosity, m 2 /s 

u turbulence frequency, 1 /s, 

frequency (=^%r - ,!>.)- 1 /s, 
complex growth rate, 1 /s, or 
acentric factor 

fl maximum growth rate, 1 /s 

p density, kg/m 3 

pin liquid density (at 1 bar, 273.15 K), kg/m 3 
p r reduced density ( =p/ p c ) 

a potential length constant, Angstrom (=0.1 nm), 
surface tension, kg/s 2 , or characteristic 
diameter of a molecule in Table 1 
r stress tensor term, kg/ms 2 , or 

characteristic breakup time (=1 /Kb), s 

6 void fraction, or 

spray cone angle, deg. 

Subscripts 

A A-th species component 
b breakup conditions 

B B-th species component 
c critical conditions 

/ fuel 

g gas or global 

i i-th coordinate, i-th species 

component, summation index, or 
imaginary component 
inj injector 

j j-th coordinate, j-th species 

component, or summation index 
k k-th droplet group, k-th coordinate, 
summation index, or liquid conditions 
associated with k-th droplet group 
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I liquid 

L liquid ligament 

in multi-component mixture 

n normal, or nth-face of the 

computational cell 
o initial conditions, 

orifice exit conditions, 
injector initial conditions, or 
oxidizer 

p particle, drop, or conditions 

associated with a grid cell 
r reduced, radial coordinate, real 

component, or a reference value 
RT Rayleigh-Taylor 

s droplet surface, or adjacent 

computational cell 
t time 

x axial or x-coordinate 

y y-coordinate 

z z-coordinate 

, partial differentiation with respect 

to the variable followed by it 

Superscripts 

mean, or average 

first order differentiation, or 

flow rate 

second order differentiation 
// fluctuations 
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1 INTRODUCTION 

There are many occurrences of sprays in a vari- 
ety of industrial and power applications and materi- 
als processing [1], A liquid spray is a two phase flow 
with the gas as the continuous phase and the liquid 
as the dispersed phase in the form of droplets or lig- 
aments [1]. The coupling between the two phases oc- 
curs through the exchanges of mass, momentum, and 
energy involving a wide range of thermal, mass, and 
fluid dynamic factors. A number of finite-difference 
formulations have been advanced over the years for 
predicting the flow (mass and momentum) and ther- 
mal properties of a rapidly vaporizing spray. Some 
of the pros and cons of various formulations can be 
found in [1-5]. 

The state of the art in multi-dimensional com- 
bustor modeling, as evidenced by the level of sophis- 
tication employed in terms of modeling and numeri- 
cal accuracy considerations, is also dictated by the 
available computer memory and turnaround times 
afforded by present-day computers. With the aim 
of advancing the current multi-dimensional com- 
putational tools used in the design of advanced 
technology combustors, a solution procedure is de- 
veloped that combines the novelty of the coupled 
CFD/spray/scalar Monte Carlo PDF (Probability 
Density Function) computations on unstructured 
grids with the ability to run on parallel architec- 
tures. In this approach, the mean gas-phase velocity 
and turbulence fields are determined from the solu- 
tion of a conventional CFD method, the scalar fields 
of species and enthalpy from a modeled PDF trans- 
port equation using a Monte Carlo method, and a 
Lagrangian-based dilute spray model is used for the 
liquid-phase representation. 

The gas-turbine combustor flows are often char- 
acterized by a complex interaction between various 
rate-controlling processes associated with turbulent 
transport, mixing, chemical kinetics, evaporation and 
spreading rates of spray, convective and radiative 
heat transfer, among others [10]. The phenomena 
to be modeled as controlled by these processes often 
interact with each other at various disparate time and 
length scales. In particular, turbulence plays an im- 
portant role in determining the rates of mass and heat 
transfer, chemical reactions, and liquid-phase evapo- 
ration in many practical combustion devices. The in- 
fluence of turbulence in a diffusion flame manifests it- 
self in several forms, ranging from the so-called wrin- 
kled or stretched flamelets regime to the distributed 


combustion regime, depending upon how turbulence 
interacts with various flame scales [69-70]. 

Most of the turbulence closure models for reac- 
tive flows have difficulty in treating nonlinear reaction 
rates [69-70]. The use of assumed shape PDF meth- 
ods was found to provide reasonable predictions of 
pattern factors and NO.y emissions at the combustor 
exit [71]. However, their extension to multi-scalar 
chemistry becomes quite intractable. The solution 
procedure based on the modeled joint-composition 
PDF transport equation has an advantage in that it 
treats the nonlinear reaction rates without any ap- 
proximation. This approach holds the promise of 
modeling various important combustion phenomena 
relevant to practical combustion devices such as flame 
extinction and blow-off limits, and unburnt hydrocar- 
bons (UHC), CO, and NCFy predictions [71]. 

With the aim of demonstrating the viability of 
the PDF approach to the modeling of practical com- 
bustion flows, we have undertaken the task of extend- 
ing this technique to the modeling of sprays as a part 
of the NCC (National Combustion Code) develop- 
ment program [7-15]. In order to facilitate large-scale 
combustor computations, we have extended our pre- 
vious work on the combined CFD/spray/PDF com- 
putations to parallel computing [10-15]. The use 
of parallel computing offers enormous computational 
power and memory as it can make use of hundreds 
of processors in concert to solve a complex problem. 
The trend towards parallel computing is driven by 
two major developments: the widespread use of dis- 
tributed computing and the recent advancements in 
MPPs (Massively Parallel Processors). The solver is 
designed to be massively parallel and automatically 
scales with the number of available processors. A cur- 
rent status of the use of the parallel computing in tur- 
bulent reacting flows involving sprays, scalar Monte 
Carlo PDF, and unstructured grids was described in 
Ref. [11]. It outlined several numerical techniques 
developed for overcoming some of the high computer 
CPU-time and memory-storage requirements associ- 
ated with the use of Monte Carlo solution methods. 
The parallel performance of both the PDF and CFD 
modules was found to be excellent but the results 
were mixed for the spray computations showing rea- 
sonable performance on massively parallel computers 
like Cray T3D; but its performance was poor on the 
workstation clusters [10]. In order to improve the 
parallel performance of the spray module, two differ- 
ent domain decomposition strategies were developed 
and the results from both strategies were summarized 
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[10-14], 

It is also well known that considerable effort 
usually goes into gridding up of complex gas-turbine 
combustor geometries. In order to allow repre- 
sentation of complex geometries with relative ease, 
we have extended our previous work on the com- 
bined CFD/spray/PDF computations to unstruc- 
tured meshes [11-15], The grid generation time asso- 
ciated with gridding up practical combustor geome- 
tries, which tend to be very complex in shape and 
configuration, could be reduced considerably by mak- 
ing use of existing unstructured grid generators. The 
solver accommodates the use of an unstructured mesh 
with mixed elements: triangular and/or quadrilateral 
for 2D (two-dimensional) geometries and tetrahedral 
for 3D. 

With the development of two computer mod- 
ules, LSPRAY - a Lagrangian spray solver [14] and 
EUPDF - an Eulerian Monte Carlo PDF solver [15], 
we were able to demonstrate the use of the joint 
scalar Monte Carlo PDF method in the modeling of 
complex multidimensional reacting flows (e.g., gas- 
turbine combustor flows) [10-15]. In this manual, we 
only concentrate on providing complete details of our 
spray module along with some more details on the 
high pressure modifications made to the gas-phase 
calculations. However, further details on the appli- 
cation of the joint scalar Monte Carlo PDF method 
to two-phase flows could be found elsewhere in [10- 
13,15]. 

In this manual, we summarize some important 
aspects of our spray formulation without making any 
attempt to provide an in-depth review on the fluid 
dynamic and transport behavior of reacting sprays 
[6-15]. Depending on the nature of the spray, an 
appropriate selection could be made from the choice 
of various well-known spray formulations (multicon- 
tinua, discrete-particle, or probabilistic) based on ei- 
ther a Lagrangian or an Eulerian representation for 
the liquid-phase equations by making use of appro- 
priate droplet sub-grid models. The present solu- 
tion procedure could be used within the context of 
both multicontinua and probabilistic spray formula- 
tions, as it allows for resolution on a scale greater 
than the average spacing between two neighboring 
droplets [1]. For NCC, the adopted choice for the 
gas phase was an Eulerian scheme. The liquid-phase 
equations form a system of hyperbolic equations and 
they could be solved by means of either an Eulerian or 
a Lagrangian representation. A Lagrangian scheme is 
chosen as it reduces the errors associated with numer- 


ical diffusion. The liquid-phase formulation is based 
on various well-established models for droplet drag; 
the vaporization models of a poly disperse spray take 
into account the transient effects associated with the 
droplet internal heating and the forced convection ef- 
fects associated with droplet internal circulation; and 
it employs models for gas-film valid over a wide range 
of low to intermediate droplet Reynolds numbers [7] . 
And the present stochastic particle tracking method 
is applicable for flows with a dilute spray approxima- 
tion where the droplet loading is low. The numeri- 
cal method could be used within the context of both 
steady and unsteady calculations [8-13]. Not consid- 
ered in the present release of the code are the effects 
associated with droplet/shock interaction and dense 
spray effects. 

Currently, most of the finite-difference tech- 
niques used for predicting the two-phase flows make 
use of the physics derived from single-component liq- 
uid droplet studies with constant properties. How- 
ever, it is well known that most of the gas-turbine fu- 
els are multicomponent mixtures of many compounds 
with a wide distillation curve [7,18]. The multicom- 
ponent nature of the liquid sprays is becoming evi- 
dent with the increasing need to use jet fuels derived 
from heavier petroleum compounds. The gasifica- 
tion behavior of a multicomponent fuel droplet may 
differ significantly over that of a pure single compo- 
nent fuel droplet [18]. Also, the calculation of the 
variable thermo-transport properties of the liquid- 
mixtures becomes more important at high pressures. 
The flame ignition characteristics (such as the phe- 
nomena associated with flame blow-off and extinction 
conditions) could also be influenced by the nonuni- 
form concentration of the fuels with different volatili- 
ties. However, the importance of the multicomponent 
liquid fuels with variable properties received little at- 
tention in the modeling of comprehensive gas-turbine 
combustor spray computations. With this in mind, 
we have made the following improvements: 

(a) In order to deal with the gas turbine fu- 
els that are mixtures of many compounds, we have 
extended the spray formulation to multicomponent 
liquid sprays. This implementation also takes into 
account the effect of variable liquid properties. 

(b) In order to take into account the nonideal 
gas behavior under high pressure conditions, we have 
completed the integration of the Peng-Robinson EOS 
into our CFD module. Also, we have completed the 
implementation of a high-pressure correction into the 
calculation of transport properties in the gas phase. 


NASA/CR— 2012-217294 


7 



These modifications would enable the calculation of 
high-pressure flow properties in the gas phase. 

In order to reduce some uncertainty associated 
with the specification of initial droplet conditions, we 
have undertaken the task of integrating an atomiza- 
tion module into our spray solution procedure. The 
atomization module was developed by CFDRC Inc. 
[47] in collaboration with the university of Wisconsin 
(UW) [45-46]. Our computational experience with 
the integration of the atomization module was sum- 
marized in [59-60]. A complete description of all the 
primary atomization as well as the secondary droplet 
breakup models contained in the CFDRC/UW atom- 
ization module is provided in the appendix. 

An understanding of flash atomization and va- 
porization behavior at superheat conditions is iden- 
tified to be a topic of importance in the design of 
modern supersonic engines, scramjet, and ramjet af- 
terburners [61-64, & 74]. It is because higher heat- 
load requirements may lead to the use of the same 
fuel as a coolant and also sometimes the nozzles are 
required to operate under low back pressures. Also, 
under normal gas-turbine operating conditions, it is 
estimated that a small fraction of the liquid fuel may 
be released by flash boiling [62]. There are also some 
reported cases of flash related engine performance 
problems in gasoline direct-injection internal combus- 
tion engines [62]. Because of its importance in some 
NASA-related applications, we have undertaken an 
assessment study to establish baseline accuracy of 
various existing CFD models used in the modeling 
of a superheated spray. As a part of this effort, we 
have incorporated a solution procedure based on the 
modeling approach developed by [62-64, & 74], 

Some important aspects of the spray module are 
summarized below: 

• An efficient particle tracking algorithm was de- 
veloped and implemented into the Lagrangian 
spray solver in order to facilitate particle move- 
ment in an unstructured grid of mixed elements. 

• LSPRAY-IV is currently coupled with an un- 
structured flow (CFD) solver of NCC [21-23], 
and an Eulerian-based Monte Carlo probability 
density function solver - EUPDF [15], which were 
selected to be as the integral components of the 
NCC cluster of modules. EUPDF provides the 
solution for the scalar fields of species and en- 
thalpy from a modeled PDF transport equation 
using a Monte Carlo method. 


• The spray solver receives the mean velocity and 
turbulence fields from the flow solver. The so- 
lution for the scalar (energy and species) fields 
could be provided by means of either a conven- 
tional CFD solver or a Monte Carlo PDF solver 
depending on the choice of the solver. 

• The spray solver supplies the spray source-term 
contributions arising from the exchanges of mass, 
momentum and energy with the liquid-phase 
to the flow solvers (CFD and/or Monte Carlo 
PDF). This information could be used in either 
conservative or non-conservative finite-difference 
formulations of the gas phase equations. 

The furnished code demonstrates the success- 
ful methods used for parallelization and coupling of 
the spray to the flow code. First, complete details of 
the spray solution procedure is presented along with 
several other numerical issues related to the coupling 
between the CFD, LSPRAY-IV, and EUPDF solvers. 
It is followed by a brief description of the combined 
parallel performance of the three solvers (CFD, EU- 
PDF, and LSPRAY-IV) along with a brief summary 
of the validation cases. 

2 GOVERNING EQUATIONS FOR GAS 
PHASE 

Here, we summarize the governing gas-phase conser- 
vation equations in Eulerian coordinates [1]. This 
is done for the purpose of identifying the interphase 
source terms arising from the exchanges of mass, mo- 
mentum, and energy with the liquid phase. They are 
valid for a dilute spray with a void fraction of the 
gas, 0, close to unity. The void fraction is defined as 
the ratio of the equivalent volume of gas to a given 
volume of a gas and liquid mixture. 

The conservation of the mass leads to: 

[P^c],^ 4” [P^ZVi],Xi = Smlc ~ ^ difc (1) 

k 

The source term is given as a summation over dif- 
ferent classes of droplets. Each class represents the 
average properties of a different polydisperse group of 
droplets. Here, ilk represents the number of droplets 
in a given class and represents the corresponding 
mass vaporization rate. 

For the conservation of the jth species, we have: 
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3 HIGH PRESSURE EQUATION OF 
STATE 


[pVcVj],t + [pV c Uiyj] iXi - [ pV c Dy jtXi \ tXi - pV c Wj = 


Smls 


y ej n k rhk 

k 


(2) 


where 


In order to calculate the high pressure gas be- 
havior, the Peng-Robinson EOS (Equation-Of-State) 
is employed for a multi-component mixture in the fol- 
lowing form [16-17,28]: 



0 and 6j 
j 


1 


For the species conservation, the source term con- 
tains an additional variable, ej, which is defined as 
the fractional vaporization rate for species j . 

For the momentum conservation, we have: 


“t - [pbrtuyUy ] _ X j T \p^c],xi 
[@VcTij],Xj [(1 = Smlm = 

y n k m k u ki - ^2 y Pk r\ n k u k i)t (3) 
k k 

where the shear stress in Eq. (3) is given by: 

7~ij — p\Ui’Xj + U-j-Xi ] ^ /iSjj V-i.Xj 

For the momentum conservation, the first source 
term represents the momentum associated with liquid 
fuel vapor and the second represents the momentum 
change associated with droplet drag. 

For the energy conservation, we have: 


[pVch] ,t+[pV c Uih] , Xi - [ 0V c \T Xi ] iXi - [(1-0) V c X t T Xi \ tXi 

-[0V c p\ t t = Smle = Y j n k mk (hs - h,eff ) (4) 

k 

Similarly, the energy conservation has a source term 
contribution from the added liquid fuel vapor and an 
additional contribution associated with the effective 
latent heat of vaporization which accounts for the 
heat flux (loss or gain) to the droplet interior from 
the ambient. 

The main purpose of the spray solver is to cal- 
culate the source terms arising from the exchanges 
of the mass, momentum, and energy. In the case of 
NCC, it provides the calculated source terms to both 
the CFD and Monte Carlo PDF solvers. 


P = 


RT a m 

V - b m U 2 + 2 b m V - bl 


where 


a m = T,i yiyj(aiaj) 1/2 (l - %), 
brri = yi^ii 


u — 0.07780 RT ic 

~ Pic 


dj = 


0.45724i? T- 


[l + /ia,(l-T; 1 / 2 )] 2 , 


Tir = T /T ic , 


( 5 ) 


f iui = 0.37464 + 1.54226W* - 0.26992w; 2 , 

u>i is known as the acentric factor of the molecules 
which is a measure of non-sphericity of the molecules, 
and kij is known as the binary interaction coefficient. 
However, the Peng-Robinson EOS is rewritten as a 
cubic EOS in terms of the compressibility factor, Z 
(= before it is solved: 

Z 3 - (1 - B*)Z 2 + (A* - 2B* - 3 B* 2 )Z 

- A*B * + B* 2 + B* 3 = 0 (6) 

where 

A* = and B* = 

We have chosen the Peng-Robinson EOS be- 
cause of its simplicity and, more importantly, it 
proved to be very useful in the supercritical droplet 
vaporization studies of [19-20]. 

Table 1 lists various physical constants for some 
of the species that we found to be of interest in our 
spray computations. It contains values for the boiling 
temperature at normal pressure, the critical tempera- 
ture, the critical pressure, the critical density, the la- 
tent heat of vaporization at normal pressure, the crit- 
ical molar volume, the molar volume at normal pres- 
sure, and the acentric factor of the molecules. Most 
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Table 1. Physical constants. 

Species 

T n b 

T c 

Pc 

Pc 

^kn 

K c 

v n 

LO 

a 

a / k 


(K) 

(K) 

(atm) 

(Kg/m 3 ) 

(KJ/kg) 

(cm 3 /g-mole) 

(cm 3 /g-mole) 


(A °) 

(K) 

CqHi4 

341.9 

507.4 

30.0 

660.0 

334.8 

370.0 

140.06 

0.296 

5.949 

399.3 

C-jHiq 

371.6 

540.2 

27.0 

682.0 

316.3 

432.0 

162.00 

0.351 

6.297 

419.031 

CsHis 

398.82 

568.8 

24.6 

718.5 

301.3 

492.0 

188.8 

0.394 

6.62 

488.15 

C10H22 

447.3 

617.6 

20.8 

728.3 

276.1 

603.0 

233.68 

0.490 

7.16 

540.06 

C\ 2 H 2 q 

489.5 

658.3 

18.0 

748.0 

256.3 

713.0 

278.54 

0.562 

7.655 

583.68 

C14H30 

526.7 

694.0 

16.0 

763.0 

240.1 

830.0 

326.62 

0.679 

8.067 

629.08 

n 2 

77.4 

126.2 

33.9 

807.1 

197.6 

90.1 

31.87 

0.039 

3.681 

91.5 

0 2 

90.2 

154.6 

50.4 

1135.7 

212.3 

74.4 

26.08 

0.025 

3.433 

113.0 

C 0 2 

00.0 

304.1 

73.8 

000.0 

000.0 

94.0 

33.32 

0.239 

3.996 

190.0 

h 2 o 

373.2 

647.3 

221.2 

958.1 

2257.2 

57.1 

19.76 

0.344 

2.641 

809.1 


Table 2. Binary Interaction Parameters, 


C 7 H 1 Q 

0=1) 

o 2 

0=2) 

n 2 

0=3) 

c o 2 
0=4) 

h 2 o 

0=5) 

C 7 H \ 6 

(i=l) 

0.0000 

0.1321 

0.1440 

0.1000 

0.1484 

o 2 

(i=2) 

0.1321 

0.0000 

-0.0119 

-0.0289 

0.0910 

n 2 

(i=3) 

0.1440 

-0.0119 

0.0000 

-0.0170 

0.1030 

C0 2 

(i=4) 

0.1000 

-0.0289 

-0.0170 

0.0000 

0.1200 

h 2 o 

(i=5) 

0.1484 

0.0910 

0.1030 

0.1200 

0.0000 
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of this data is collected from [16-18] and is useful in 
both the evaluation of the Peng-Robinson EOS and 
the calculation of various other variable properties. 
And Table 2 lists the binary interaction parameters, 
kij , used in a calculation for the multi-component 
mixture of n-heptane, O 2 , N 2 , CO2, & H 2 0. While 
the data for most of the binary pairs was obtained 
from various reference books, some of the missing 
data, however, is replaced with known data found for 
other binary pairs of molecules with similar molecular 
weights. 

4 HIGH PRESSURE CORRECTIONS FOR 
GAS-PHASE TRANSPORT PROPERTIES 

The effect of pressure on the viscosity of pure 
gases is determined by making use of the Reichenberg 
method [16,19,29]: 



0.0 1.0 2.0 3.0 4.0 5.0 6.0 

Reduced pressure, P r 


Figure 1 Takahashi correlation showing the 
variation of the binary diffusion coefficient 
versus reduced pressure at different reduced 
temperatures. 


Pn 


1 + Q v 


A V P, 


, 3/2 


B v P r + (1 + C v P r U ”) 


*M-1 


( 7 ) 


where 


A v = exp(a 2 T . ’“) B v = A v {j3{T r - fl 2 ) 

C v = ^exp(7i T r c ) D v = j^exp(\ 2 Tr) 


ol\ = 1.982510-' 
/3i = 1.6552 
71 = 0.1319 
Ai = 2.9496 


a 2 = 5.2683 
P 2 = 1.2760 
72 = 3.7035 
A 2 = 2.9190 


a = -0.5767 

c = -79.8678 
d= -16.6169 


and Q v = 1.0 for non-polar molecules. 

The effect of pressure on the thermal conduc- 
tivity of pure gases is determined by making use of 
the Stiel & Thodos method [16,19,30]: 


(A - A n )r Z\ = 1.2210 _2 [exp(0.535/9 r ) - 1] (8) 

when p r < 0.5 

(A - A„)rZ c 5 = 1.1410 _2 [exp(0.67 / o r ) - 1.069] (9) 
when 0.5 < p r < 2.0 

(A - A„)rZ c 5 = 2.6010- 2 [exp(1.155p r ) + 2.016] (10) 
when 2.0 < p r < 2.8 


where A is in W/(m.K), T = 210 [ Tl pf ] 1 ^ 6 , Z c is the 
critical compressibility, and p r is the reduced density 

p/Pc = Vjv. 

The effect of pressure & temperature on diffu- 
sion coefficients in gases is determined by making use 
of Takahashi correlation [16,19,31]: 


_D_abP_ 
(■ D ab P) + 


f(T r ,Pr) 


( 11 ) 


where Dab = diffusion coefficient, cm 2 /s, P = 
pressure, bar, the superscript + indicates that low- 
pressure values are to be used, and the reduced pres- 
sure and temperature in the above equation are cal- 
culated as follows: 


T r 

T c ab 

P r 

PcAB 


T 

TcAB 


UaTcA + UbTcB 
p 


PcAB 

UaPcA + VbPcB 


This correlation is valid for a system involving 
a trace solute diffusing in a supercritical fluid. It is 
shown to provide satisfactory results but it is based 
on a very limited set of available experimental data. 
The graphical representation of the Takahashi func- 
tion is given in Fig. 1. 

It is noteworthy that at low pressures, the poly- 
nomial fits for the variable thermodynamic properties 
are taken from the data set compiled by McBride et el 
[24]. The transport properties involving the thermal 
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conductivity and molecular viscosity for individual 
species is estimated based on the Chapman-Enskog 
collision theory [25-27] . And Wilke’s formulae is used 
to determine the properties of mixture [25-27]. The 
binary-diffusion coefficients are determined based on 
the Chapman-Enskog theory and the Lennard-Jones 
potential [25-27]. 

5 LIQUID-PHASE EQUATIONS 

Here, we summarize the governing equations for 
the liquid-phase based on a Lagrangian formulation 
where the equations for particle position and velocity 
are described by a set of ordinary differential equa- 
tions. For the particle position of the kth droplet 
group, we have: 


dx 


ik 


dt 

For the droplet velocity: 


— ^ ik 


(12) 


dU'ik 3 C'D^gsRe-k r , / 1 /i 

~dT= 16 ^2 K + u g ~ u ‘k\ (13) 


p k r% 


where 


2 / 3 ' 


c — 24 | i I 

Co ~lte k [ 1+ -6- 


(14) 


According to Yuen and Chen [32] , the droplet drag for 
a vaporizing droplet could be calculated by a solid- 
sphere drag correlation but suggested using a correc- 
tion in the evaluation of the droplet Reynolds number 
for the average viscosity based on a 1/3 weighting rule 
as given by Eq. (46). The droplet Reynolds number 
is given by: 


Re k = 2 7 - A ^ [( u g + u' g - u k ) . ( u g + u' g - u k )] 1/2 

Pgs 

(15) 

By following the approach taken from KIVA 
[33], a gas turbulence velocity, u' g , is added to the 
local mean gas velocity when calculating the droplet 
drag and vaporization rates. The gas velocity fluctu- 
ations is calculated by randomly sampling a Gaussian 
distribution with mean square deviation, 2/3 k. The 
Gaussian is given by 


G(u' g ) = (4/37rfc) 3 / 2 exp[—3\u' g \ 2 /4k\ (16) 


The gas fluctuating component is calculated once at 
every turbulence interaction time, t tur , and is oth- 
erwise held constant [33]. The correlation time is 
taken to be the minimum of either the eddy time or 
the transit time taken by the droplet to traverse the 
eddy. It is given by 


. (k k 3 1 2 1 \ 

ttur — TTl'lTl ( , Ctt "j ; r I 

Ve e \u g + u'g - u k \J 


(17) 


where Ctt is an empirical constant with a value of 
0.16432. 

The liquid mixture density, p k , in Eq. (13) is 
given by [16,18]: 


Pk — ^ ) VkiPki (13) 

i 

and the individual component liquid density is given 

by: 


M ki 
Pki ~ Vm 

where the molar volume, V k i, is 

V ki = U C i(0. 29056 - 0.8775wi) c ’ 

and 


(19) 


(20) 


The droplet regression rate is determined from 
one of three different correlations depending upon the 
droplet Reynolds number range. The first correlation 
as given by Eq. (21) is based on a gas-film analysis 
developed by Tong & Sirignano [34]. It is based on a 
a combination of stagnation and flat-plate boundary- 
layer analysis and is valid for Reynolds numbers in 
the intermediate range. The last two correlations as 
given by Eqs. (22) & (23) are taken from Clift et 
al [35] and are valid when Re k < 20. In fact, when 
the droplet Reynolds number goes to zero, Eq. (23) 
becomes identical to the droplet regression rate of a 
vaporizing droplet in a quiescent medium [27] . 


ds k 

dt 



r 2 I 1 / 2 
—Re k f{B k ) 

7 r 


if Re k > 20 


(21) 
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ds k 

dt 


PgDg 


Pk 

if 1 < Re k < 20 


1 + (1 + Rek) 1 / 3 


Rel 077 


ln( 1 + B k ) 

(22) 


the droplet regression rate. The initial and boundary 
conditions for Eq. (27) are given by: 


t = t. 


a = 0, 


injection ? 

dTk = 1_ 
da 17 


Tk = 

Cpipi 

A; 


■ ar fc 
: at 


(29) 

(30) 


dsfc 

dt 


Pg D g 


Pk 

if Re k < 1 


1 + (l + i?e fc ) 1/3 l in(l + B fc ) (23) 


where B k is the Spalding mass transfer number and 
is given by: 


Bk = 


(Vfs - Vf) 
(! - Vfs) 


(24) 


and yf s is given as a summation over the fuel-species 
mass fractions at the droplet interface: 


yfs = E yi 


(25) 


and yf is given as a summation over the fuel-species 
mass fractions in the ambient: 


Vf 


E 


Vfi 


(26) 


9Tk_ A , 

dt Cpipirl |P da 2 


d 2 T fe \ dTk 

a - 0 + (1 + C(t)a) — 


5a 


where 


CM = T7 


C plP i 
A i 


r k 


dr k 

dt 


dT k 3 p k ds k 

a = 1) a^ = “32A; [ ^ e// “ Zfcl lt- (31) 

where a = 0 refers to the vortex center, and a = 1 
refers to the droplet surface, and the mixture latent 
heat of vaporization, l k , is given by 


ifc — ^ ^ 


(32) 


and the individual component latent heat of vapor- 
ization, l k i, is given by [16,18]: 


Iki — 


Ikin ^ 


T c i ~ T k 


0.38 


(33) 


and the function f(B k ) is similar to that of a Blassius 
function [1, 36] and is obtained from an analysis simi- 
lar to that of Emmon’s boundary-layer flow over a flat 
plate with blowing [36]. The range of validity of this 
function was extended in Raju and Sirignano [7] to 
consider the effects associated with a boundary-layer 
flow with suction. 

The internal droplet temperature is determined 
based on a vortex model [34] . The governing equation 
for the internal droplet temperature is given by: 


" T ci — T bi 

and the droplet boiling temperature is given by 

rp IkinMi/ R u 

bl lkinM i /(R u t bni ) - ln(P) 

and, finally, the effective latent heat of vaporization, 
lk,efft is defined as: 


(34) 


lk,eff = lk+ 47^ 

m. k V dr J , 


(35) 


It is a very useful parameter as it represents the total 
energy loss associated with the latent heat of vapor- 
ization in addition to the the heat loss to the droplet 
interior. l k ,eff is calculated by means of the following 
relationship [18]: 


lk,eff — 


C p (T g ~T ks ) 

(i + Bk y/^-i 


(36) 


Similar to the internal thermal transport, the 
internal mass transport of a multi-component fuel is 


(27) 

given by: 

(28) 

dy ki 

dt 


,D k 

77T 


' k L 


d 2 Vki . n , dy ki 

a - „ + (1 + 6 (t)a) 


da 2 


da 


(37) 


where a represents the coordinate normal to the 
stream-surface of a Hill’s Vortex in the circulating 
fluid, and C{t) represents a nondimensional form of 


The initial and boundary conditions for Eq. (37) are 
given by: 
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t = t 

injection-) 

Dk 

i Vki,o 

(38) 

a = 0, 

dyki _ 

1 

r r 2 1 

' k 

dyki 

(39) 


da 

17 

_D k _ 

dt 

, dy ki 

a = l, ~ 

aa 

32 D k ^ Vkis 

, ds k 

(40) 


By knowing the mass fractions of the liquid 
species at the droplet surface, the corresponding mole 
fractions are determined by 


%iks 


Viks / 'Mi 

X)i 2 lyka/Mi 


(41) 


At the droplet interface, the mole fractions of 
the gas species are obtained by means of Raoult’s 
law: 


%i s p^iks^is (42) 

where the partial pressure, P, s , is determined by 
means of the Clausius-Clapeyron relationship: 


P is = exp 


Iki j 

( 1 

AX 

Ru 

{T bl 

Tk s J _ 


(43) 


Then, the corresponding mass fractions in the gas- 
phase at the droplet interface are given by 


v . = XiaMi (44) 

Uls M a (l ~ Ei x is) + Ei MiXis y ’ 

where M a is the molecular weight of the gas excluding 
fuel vapor. And the fractional mass vaporization rate 
of liquid species, ej, is given by 

£i = Vis + (1 ^ Vfs) — — — (45) 

y/s - vs 

It is noteworthy that the thermodynamic and 
transport properties at the gas film are calculated at 
the temperature and composition as determined by 
the following one-third rule: 

1 2 

( fiavg — (46) 

The correlations for the gas-phase thermody- 
namic and transport properties are described in Sec- 
tion 4. In a similar way, the liquid-phase thermody- 
namic and transport properties are determined based 
on the correlations described in the next section. 


6 LOW PRESSURE LIQUID 
THERMODYNAMIC & TRANSPORT 
PROPERTIES 

The specific heat at constant pressure, C p i, 
thermal conductivity, A/, and viscosity, f. q, are evalu- 
ated by means of the following expressions: 


Cpi — C p io + CpnT + CpiiT" + CpizT 3 + Cpi^T 4 ^ 

(47) 


A; = A,o + X n T + A l2 T 2 + A l3 T 3 + \ U T 4 + A , 5 T 5 

(48) 


In pi = pw + pn/T + p l2 T + pi 3 T 2 (49) 


where C p i is in J/(kg K), pi in (pPA s), and A i in 
(W/m K). 

Tables 3-5 provide the polynomial constants 
used in Eqs. (47)-(49) for some of the species that 
we found to be of interest in our spray computations. 
Table 3 provides the constants for the liquid specific 
heat, C p i, Table 4 for the liquid thermal conductivity, 
A i, and Table 5 for the liquid molecular viscosity, /q. 
These tables are compiled with the data taken mostly 
from the references of [16,18]. 

The binary diffusion coefficient, Dij, is evalu- 
ated as follows [16]: 


— 


KdifT 


(50) 


where Kdif = 8.210 _s (l + [qt 1 ] 2 ^ 3 ). One should be 
careful in using this approximation as it is based on 
a scarce set of available experimental data. 

The specific heat for a multicomponent mixture 
is given by: 


n 

c pm = J2y' c p* ( 51 ) 

i= 1 

And the thermal conductivity of a multicomponent 
mixture is calculated by means of the Li method [16]. 

n n 

A m = 'y ' y ' ( 62 ) 

*= 1 3 = 1 

where 
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Table 3. Polynomial constants for liquid specific heat. 

Fuel 

CpZO 

CpZl 

Cpl2 

Cpl 3 

CpZ 4 

C(>H\4 

2.4169 

-5.9866e-03 

2.0959e-05 

-8.4489e-09 

0.0 

C7H1Q 

4.8227 

-3.6980e-02 

1.6777e-04 

-3.0987e-07 

2.2081e-10 

CgHig 

9.2189 

-8.8314e-02 

3.8869e-04 

-7.2539e-07 

5.0776e-10 

C 10 H 22 

4.7991 

-2.8643e-02 

9.3619e-05 

-8.9516e-08 

0.0 

C 12 H 26 

4.7900 

-2.8643e-02 

9.3619e-05 

-8.9516e-08 

0.0 

Cl4^30 

4.7991 

-2.8643e-02 

9.3619e-05 

-8.9516e-08 

0.0 



Table 4 

. Polynomial constants for liquid thermal conductivity. 


Fuel 

Azo 

Azi 

Az2 

A/3 

A/4 

A/5 

CqH\4 

0.37078 

-5.4313e-03 

4.628e-05 

-1.8002e-07 

3.2243e-10 

-2.1832e-13 

C7H1Q 

0.13236 

9.4441e-04 

-6.588d-06 

1.4617e-08 

-1.1244e-ll 

0.0 

CsHis 

0.25652 

-7.5401e-04 

1.5872e-06 

-1.6795e-09 

-1.3375e-16 

0.0 

C 10 H 22 

0.22179 

-2.3699e-04 

-6.94e-07 

2.0415e-09 

-1.5741e-12 

0.0 

C 12 H 26 

0.17609 

4.2463e-05 

-7.4467e-07 

6.9446e-10 

0.0 

0.0 

C 14 H 30 

0.18801 

-9.1399e-05 

-2.1464e-07 

1.1655e-10 

0.0 

0.0 

n 2 

-2.629E-1 

-1.545E-3 

-9.450E-7 

0.0 

0.0 

0.0 

0 2 

2.444E-1 

-8.813E-4 

-2.023E-6 

0.0 

0.0 

0.0 

C0 2 

4.070E-1 

-8.438E-4 

-9.626E-7 

0.0 

0.0 

0.0 

h 2 o 

-3.838E-1 

5.254E-3 

-6.369E-6 

0.0 

0.0 

0.0 


Table 5. 

Polynomial constants for liquid molecular viscosity. 

Fuel 



H 12 


CqH \4 

-4.034E+00 

8.354E+02 

0.0000000 

0.0000000 

C 7 H 1 Q 

-4.325E+00 

1.006E+03 

0.0000000 

0.0000000 

CgHig 

-4.333E+00 

1.091E+03 

0.0000000 

0.0000000 

C 10 H 22 

-4.460E+00 

1.286E+03 

0.0000000 

0.0000000 

C 12 H 26 

-4.562E+00 

1.454E+03 

0.0000000 

0.0000000 

C 14 H 30 

-4.615E+00 

1.588E+03 

0.0000000 

0.0000000 

n 2 

-2.795E+01 

8.660E+02 

2.763E-01 

-1.084E-03 

0 2 

-4.771E+00 

2.146E+02 

1.389E+02 

-6.255E-05 

C0 2 

-3.097E+00 

4.886E+01 

2.381E-02 

-7.840E-05 

h 2 o 

-2.471E+01 

4.209E+03 

4.527E-02 

-3.376E-05 
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Xij — 2(A, ; 1 + X J 1 ) 1 

A. — XiVi 

where x-i is the mole fraction of the species i, (pi is a 
volume fraction of the ith species, and V % is the molar 
volume of the pure fluid. 

7 VAPORIZATION MODELING OF A 
SUPERHEATED DROPLET 

Flashing phenomena refers to a process that is 
in thermodynamic non-equilibrium when a liquid is 
superheated [72-73]. The reasons for its occurrence 
are mainly two-fold [72-73]: (1) a liquid fuel can be 
heated to a higher temperature (above its saturation 
temperature) while its pressure is maintained, and 
(2) when a liquid is depressurized rapidly it can lead 
to flash injection because the thermal inertia initially 
tends to maintain its bulk internal temperature above 
the saturation temperature. Although flash evapora- 
tion is considered to be detrimental to engine perfor- 
mance under normal circumstances, it can have some 
potential benefits: it is known to produce a fine spray 
with enhanced atomization, increase effective spray 
cone angle, and decrease spray penetration [61]. 

An understanding of flash injection is of im- 
portance in some applications involving scramjet and 
ramjet afterburners because the same liquid fuel is of- 
ten used as a coolant coupled with engine conditions 
where nozzles operate at low back pressures and su- 
personic outflow [61]. The objective of our work is to 
establish a baseline accuracy for existing atomization 
and vaporization models used in the prediction of a 
superheated spray by undertaking a critical review 
of existing experimental data available in the litera- 
ture for validation. This work is funded through the 
supersonics (SUP) and subsonic fixed wing (SFW) 
project office initiatives on high altitude emissions of 
the NASA’s fundamental aeronautics program. 

In what follows we first provide some details of 
the two superheat vaporization models that we have 
implemented into NCC. It is followed by other model- 
ing considerations that need to be taken into account 
as described in Sections 7.3 & 7.4. 

7.1 Superheat Vaporization Model of Zuo, Gomes, & 
Rutland [62], and Schmehl & Steelant [63-64] 


gas. Under superheat conditions, the characteristic 
vaporization time resulting from the external heat 
transfer from the surrounding gas is of the same order 
of magnitude as that resulting from the flash evap- 
oration. The energy needed for vaporization at the 
droplet surface is partly provided by the superheat 
energy stored within the droplet and it is controlled 
by the droplet internal heat transfer. The modeling 
approach differs from the classical droplet vaporiza- 
tion models in three important ways: (1) under su- 
perheat conditions, the droplet surface mass fraction, 
Yf s , approaches unity as the droplet surface temper- 
ature is maintained at the corresponding liquid boil- 
ing temperature; (2) under superheat conditions, all 
the external heat transfer from the surrounding gas 
is made available to the vaporization process with no 
apparent increase in the droplet surface temperature; 
and (3) the flow of fuel vapor imparted by flash va- 
porization partly counterbalances the flow generated 
by external heat transfer and may significantly reduce 
the energy transferred from the surrounding gas. 

Based on the governing equations of conserva- 
tion for an isolated spherically symmetric droplet, 
Zuo et al [62] and Schmehl and Steelant [63-64] 
showed that the total evaporation rate, m can be 
calculated as 

irik = m k , flash + mk,t (53) 

where the flash boiled vaporization rate, rrik, flash, is 
given by 

rhkjiash = 47T rla s ^ k ^ (54) 

where T k is the internal droplet temperature and the 
overall heat transfer coefficient, a s (= kJ/s m 2 °K) 
is given by the Adachi correlation [65]: 

= 0.76(Tfc - T t ) 0 ' 26 (0 < T ks -T h < 5) 

a B = 0.027(7), - T b ) 233 (5 < T ks - T b < 25) (55) 

= 13.8(7), - T b )°' 39 ( T ks -T b > 25) 

The Adachi correlation is valid over a wide range of 
superheat conditions. 

The vaporization rate due to external heat 
transfer, rhfc jt , in Eq. (53) is given by 


It is based on an extension of the classical D 2 - 

theory. In the classical evaporation model, the ther- m k , t = 2 ^ in[ 1 + (1 + mk ^ lash )g f ] ( 56 ) 

‘ , , . C p 1 . 7n k , flash m k ' t 

mal energy needed for evaporation is mostly furnished 
by the external heat transfer from the surrounding 
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where the Spalding heat transfer number, B t , and the 
Nusselt number, Nu, are given by 

B t = Cp ( f 9 ~ Tks) (57) 

Nu = 2(1 + 0.3i?e 1/2 Pry 3 ) (58) 

The corresponding droplet regression rate, 
is given by 


lifetimes are on par with those vaporizing under non- 
superheat conditions. So its use is likely to be re- 
stricted for the modeling of smaller droplets produced 
after flash induced atomization. 

The solution for Eq. (60) and all the JP8 fluid 
properties needed in this calculation were obtained 
by the computer coding provided by UTRC. 

7.3 Vaporization Model Valid Under Non-Superheat 
Conditions 


ds k m k 

dt 27T r k pi 


(59) 


This model is valid over an entire range of su- 
perheat conditions as long as there is some amount of 
superheat energy available within the droplet (T k > 

T b ). 


7.2 Superheat Vaporization model of Lee et al. [74] 

This modeling approach is also based on an ex- 
tension of the classical _D 2 -theory but it differs from 
the previous model of [62, 63-64] in several important 
ways: (1) neglects the effects of internal nucleation 
and convection, (2) instead of invoking the Adachi 
correlation for the flash induced vaporization rate, 
the heat loss from the droplet interior is modeled by 
<7 / = T^,o^t 3 (Poc) ( 3) the formulation is valid 

only when T, ' x > Tb u b(Poo), which means for posi- 
tive Spalding heat transfer numbers. 

With the assumptions employed, Lee et al. [74] 
arrived at the following equation for the evaporation 
rate, m v . 


Under moderate initial superheat conditions, 
only a fraction of the vaporization takes place under 
superheat conditions {T k > Tb) and the remainder 
takes place under non-superheat evaporating condi- 
tions (T k < Tb). So there is a need to revert back to 
a vaporization model valid under stable evaporating 
conditions when the internal droplet temperature ap- 
proaches the boiling temperature. In the present cal- 
culations, the vaporization rate under normal evapo- 
rating conditions is evaluated by means of a simplified 
classical Z? 2 -theory: 


TO fe = 2in'kPg Df gs Sh ln( 1 + B k ) (61) 

where the Spalding mass transfer number, B k , is 
given by Eq. (24) and the Sherwood number, Sh, 
is given by 

Sh = 2(1 + 0.3 Re 1 ' 2 Sc 1 / 3 ) (62) 

7.4 Internal Droplet Temperature Calculation 


— ~T = rhyL-qi (60) 

e m ” — 1 

where 

T = C pg T / L r , m v = (Anr^pgu) / (An A g rk,t.=o / C pg ) , 
L = L k lL r . The solution for the above equation is 
obtained by employing the standard Newton’s itera- 
tion. 

It is noteworthy that because of the assump- 
tions invoked in this model, the applicability of this 
model may be limited in several important ways: 
(1) our experience shows that the assumption of 
Too > Tbub(Poo) becomes too restrictive in most 
non-reactive calculations and also in some regions of 
most reactive spray calculations, and (2) the qi term 
is modeled based on the assumption that the overall 
droplet life is very short (sub millisecond range) but 
in most moderate superheat conditions the droplet 


Our experience with the validation studies 
showed us that there is a definite need to include 
a calculation involving the internal droplet tempera- 
ture valid under both superheat and normal evapo- 
rating conditions. In our present calculations, it was 
evaluated by means of a simple infinite conductivity 
model. 


dT k _ 3 [h,eff ~ h] ds k 
dt 2C p ir"l dt 

if T k < T b , and 


(63) 


dT k 

dt 


3cr s 

T kPlT’pl 


(T k 


Tb) 


( for the superheat model of [62, 63 — 64]) 


or 
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dT k _ 3A i dT 
dt r k piC p i dr 


(64) 


( for the superheat model of [74]) 


if T k > T b 


8 DETAILS OF ATOMIZATION 
MODELING 

The atomization refers to a process of the liq- 
uid jet breakup into droplets. There are many pro- 
cesses associated with the liquid jet breakup. They 
can broadly be classified into two regimes: (1) the in- 
ner nozzle flow, and (2) the external nozzle flow. In 
the inner nozzle flow, several factors (such as injector 
type, geometry, and size) influence the conditions at 
the injector exit (such as the velocity, the initial sheet 
or jet thickness, and the angle of droplet dispersion). 
Right now, we rely on widely-used correlations in de- 
termining the initial spray conditions at the injector 
exit. For a better description, a more accurate anal- 
ysis is needed which takes into account the physics 
associated with inside bubble growth, cavity forma- 
tion, and internal turbulence. In the external noz- 
zle flow, the liquid becomes unstable under the influ- 
ence of aerodynamic instabilities and finally breakup 
into droplets. The widely known aerodynamic insta- 
bilities are of Rayleigh-Taylor and Kelvin-Helmholtz 
kind [41, 43-44]. The Rayleigh-Taylor instability is 
due to inertia of the denser fluid opposing the system 
acceleration in a direction perpendicular to the in- 
terface of the denser fluid and the Kelvin-Helmholtz 
instability is caused by the viscous forces due to the 
relative motion of the fluids [43-44] . When the max- 
imum amplitude of the most unstable wave reaches 
a critical value, some liquid is stripped of the main 
liquid core in the form of primary droplets. These 
droplets may further breakup into smaller droplets 
due to a process known as the secondary droplet 
breakup mechanism. 

Based on a linear instability analysis of a 2D 
viscous incompressible fluid moving thorough an in- 
viscid incompressible gas, Reitz and Bracco [41] char- 
acterized the break-up regimes to be four-fold: (1) 
Rayleigh breakup, (2) first wind-induced breakup, 
(3) second wind-induced breakup, (4) atomization. 
In the first two regimes, drops of sizes greater than 
or equal to the nozzle diameter are produced at dis- 
tances far from the nozzle exit. In the applications of 


our interest, the last two regimes are more important 
where droplets much smaller than the nozzle diame- 
ter are produced near the nozzle exit. The knowledge 
gained from the instability analysis of various kinds 
[41, 45-46] is combined with some experimental ob- 
servations to form the basis for both primary atom- 
ization as well as secondary droplet breakup model- 
ing of liquid sprays. In this aproach, the jet breakup 
is modeled first by making use of a drop represen- 
tation approach in which discrete parcels of liquid 
are injected in the form of blobs with a characteris- 
tic size representative of the nozzle diameter instead 
of tracking an actual intact liquid core that forms at 
the nozzle exit. In the case of a planar or conical 
liquid sheet, the discrete parcels essentially represent 
liquid ligaments. Before the jet breakup the discrete 
parcels stay inside of the liquid core or sheet but af- 
ter the jet breakup they move independently. The 
breakup criterion, atomization rate, drop size and ve- 
locity and the location of the newly formed droplets 
are primarily determined based on an instability anal- 
ysis derived from the conservation equations of mass, 
momentum and energy. The analysis of the jet or 
sheet breakup into ligaments or droplets, the strip- 
ping of the liquid into fragments or droplets, and the 
formation of more droplets from further breakup of 
ligaments or fragments are all described under the 
classification of primary atomization. 

Some of the large droplets that are formed im- 
mediately after the primary liquid jet breakup may 
further breakup into smaller droplets under the influ- 
ence of aerodynamic instabilities. The large droplets 
first tend to flatten under the influence of aerody- 
namic pressure. Then large amplitude long wave- 
length waves caused by drop deceleration induce 
a Rayleigh-Taylor instability on the flattened drop 
causing it to breakup further into several relatively 
large-size product droplets. While at the same time 
short surface waves induce a Kelvin-Helmholtz insta- 
bility on the windward side of the parent drop re- 
sulting in the generation of much smaller product 
droplets. The breakup of the larger droplets into 
smaller droplets is described under the classification 
of secondary droplet breakup. 

Most of the coding required for both the pri- 
mary atomization and secondary droplet breakup 
models was provided by CFDR.C. These models are 
applicable under both non-superheat and superheat 
conditions. In the superheated regime, the thermo- 
physical liquid properties are replaced by the two- 
phase properties of a superheated fluid [74]. The CF- 
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DRC atomization module contains the following four 
primary atomization models: (1) the sheet breakup 
primary atomization model, (2) the blob jet pri- 
mary atomization model, (3) the air-blast atomiza- 
tion model & (4) the BLS (Boundary-Layer Strip- 
ping) primary atomization model; and the follow- 
ing three secondary droplet breakup models: (1) the 
Rayleigh-Taylor secondary droplet breakup model, 
(2) the TAB (Taylor Analogy Breakup) secondary 
droplet breakup model, and (3) the ETAB (Enhanced 
TAB) secondary droplet breakup model. The choice 
between various models depends on the specific ap- 
plication. Complete details of the models contained 
in the atomization module can be found in the ap- 
pendix. 

8.1 Flash Atomization Modeling 

There are three components to the atomization 
modeling of a superheat spray: (1) primary atomiza- 
tion, (2) secondary droplet breakup, and (3) flash in- 
duced atomization. Further details of both primary 
atomization and secondary droplet breakup models 
are described in a previous section. For the super- 
heated sprays, we have to deal with a third compo- 
nent known as flash induced atomization. We have 
followed the modeling approach of Lee et al. [74] for 
the flash induced atomization. In what follows we 
describe some details of their modeling approach. 

Flash atomization is a result of rapid evapora- 
tion associated with strong nucleation and bubble 
growth occurring inside of a non-equilibrium two- 
phase fluid [74]. When the void fraction within the 
rapidly evaporating fluid reaches a certain critical 
void fraction, the inner liquid core would shatter into 
pieces. It seemed to have a value of 0.61 based on 
the research identified in [74]. Incidentally, this con- 
dition for criticality also coincides with the packing 
limit of spheres beyond which it can not accommo- 
date any more bubbles and results in shattering [75]. 
In order to calculate this state of critical limit, they 
employed the homogeneous relaxation model (HRM) 
for calculating the mass fraction of gas, x, also known 
as quality. 


e = 


Pi~ P 


Pi,(h) - P 


( 66 ) 


Pi~ Pv P c - Pb{h ) 

Once the critical limit for the void fraction is 
reached, the resulting droplets produced as a result 
of flash atomization are obtained by a droplet number 
density correlation developed by Kawano et al. [77] . 


n = Aem-T<,u b ) e at+b (67) 

From Eq. (67), the mean droplet size, R, is 
calculated by equating total liquid mass before and 
after the droplets are formed. 


i? 4 = 


1 - e c 


(68) 


(4/3)7rn 

The solution for Eqs. (65) to (68) and all the 
JP8 fluid properties needed in this calculation were 
obtained by the computer coding provided by UTR.C. 
Also, it is noteworthy that the properties for JP8 were 
predicted by means of a surrogate model of Violi. 
Further details can be found in [74], 

The flash induced atomization is combined with 
the primary and secondary droplet breakup models 
following the approach of Lee et al. [74]. The aero- 
dynamic instability mode will atomize the liquid core 
when the surface disturbance reaches a certain critical 
value and the flash induced atomization is responsi- 
ble for shattering the liquid core when the void frac- 
tion reaches the critical value. Therefore, the control 
mechanism depends on whichever mechanism is first 
responsible for reaching the critical condition. 

8.2 Some Effects of Superheat on Primary 
Atomization 


Flash evaporation is known to produce a fine 
spray with enhanced atomization, increase effective 
spray cone angle, and decrease spray penetration [61]. 
Here, we consider the effects of flash evaporation on 
the sheet breakup primary atomization model by fol- 
lowing the approach of Zuo et al [62]. Its effect on 
the initial droplet size generated immediately after 
the first ligament breakup, d )g , is given as a function 
of both engine pressure and a superheat parameter. 


dx 

dt 


0 = e n 


(65) 


= di n ( 


Pat 


-) [i-x(- 


] o < X (- 


’ < 1 
( 69 ) 


where the values for a (= —0.54), (3 (= —1.76), and 
0 o (= 6.51.10 04 s) are taken from [76], and the void 
fraction, e, and the degree of superheat, </>, are defined 
as follows: 


where di n is the corresponding droplet droplet size 
occurring under normal evaporating conditions with- 
out flash evaporation, and y, a superheat parameter, 
is defined as follows: 
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I(T k ) - I(T b ) 
l(T b ) 


(70) 


initial velocity variation within different droplet 
classes arising from the geometric considerations 
of a chosen spray, & 


where / is the internal energy of the liquid. Its value 
varies between 0 < x < 1 with \ = 0 referring to zero 
flash evaporation and x = 1 to full flash evaporation. 

In Eq. (69), the increase in d, s due to an in- 
crease in engine pressure by a factor of (py— ) 0 ' 27 
is based on an experimental correlation obtained 
from Lefebvre [66]. It reflects the influence of cham- 
ber pressure on wave propagation as it damps wave 
growth. But the decrease by a factor of (1 — 
X(5pm)°- 135) j s due a significant reduction in 
droplet size caused by both cavitation and bubble 
growth under flash evaporation conditions. It was 
introduced based on the experimental data obtained 
from VanDerWege et al [67] and Reitz [68]. 

As the liquid approaches boiling, it also causes a 
substantial decrease in both intact liquid core length 
and core droplet size leading to a modification of the 
nominal cone angle, 0, as given by 

e = o n + (144 -o n )x 2 (71) 

where 9 n is in degrees for a spray vaporizing under 
normal conditions without flash evaporation. This 
correlation was developed based on the experimental 
data of Reitz [68] . 

These two correlations were developed in con- 
junction with the sheet breakup primary atomization 
model but their application with other primary at- 
omization models need further investigation. 

9 DESCRIPTION OF INITIAL SPRAY 
CONDITIONS 

The spray computations facilitate the use of 
multiple fuel injectors. The same or a different type 
of liquid fuel can be specified for each one of different 
injectors. The initial droplet temperature is assumed 
to be the same for all different droplet groups of any 
given injector. The liquid fuel injection is simulated 
by introducing a number of discretized parcels of liq- 
uid mass at the beginning of every fuel-injection time 
step, A tu- The following three variables play an im- 
portant role in simulating the injector initial condi- 
tions: 

• no_of_holes() - the number of holes per injector, 

• no_of_streams() - the number of droplet streams 
per hole - it is introduced to distinguish the 


• no_of_clroplet_groups() - the number of droplet 
groups per stream - it is introduced to distin- 
guish the droplet-size variation within different 
droplet classes of a polydisperse spray. 

The total number of droplet groups introduced at 
the beginning of every different injection time step 
for a given injector is therefore equal to a value 
given by the multiplication of these three vari- 
ables. Further details on some of the spray in- 
put variables can be found in the spray input file, 
ncc Jnjector.in.l of Table 6, where the integer num- 
ber followed by the last dot represents the in- 
jector number, which in this case happens to be 
one. The table provides a complete description of 
the following input variables: nolc(n_i), out_string, 
(ymki(n j,n_l), nd=l,nolc(nJ)), tdrop(n_i), atomiza- 
tion(n_i), drop .breakup _model(nJ), spray _table(n J), 
steady spray model. no_of_holes(n J) , 

no_of_streams(nd), 

no_of_droplet_groups(n j), lmdis(n_i), smdm(nj), 
cone(nJ), size_min(n_i), size_max(nJ), stochas- 
tic(n_i), ((x_inj(nJ,nx), y_inj(nJ,nx), z jnj(n_i,nx), 
z_inj(nJ,nx), flowf(n_i,nx), vJnj(n_i,nx), 

alpha inj(n i,nx), betadnj(n_i,nx), theta inj(n i,nx), 
dtheta_inj (n _i,nx) , swlr .angle (n_i,nx) ) , 

nx= 1 , no _of .holes (n J) ) , ( at om_ty pe (n _i , nx) , 

breakup _type(nd,nx), dia_hole(n_i,nx), 

delp inj(n i,nx), liq_vel(n j,nx), gas_u(n _i,nx), 

gas_v(nJ,nx), gas_w(n _i,nx), pcl_start(n _i,nx), 

pcLend (n J ,nx) ) , nx= 1 ,no_of -holes (n_i) ) . 

The initial droplet distribution for a given in- 
jector could be specified by making use of one of the 
three available options: (1) by providing a complete 
specification of the initial conditions through a spray 
table, (2) by invoking certain pre-defined correlations, 
or (3) by chosing one of four available primary at- 
omization models. When Option (3) is selected, the 
logical variable, atomization(), is set to .true, or oth- 
erwise it is set to .false.. 

Option (1): 

When a spray table is defined, one should pro- 
vide the information on the (x,y z) components of 
the initial droplet location, the (u, v, w) compo- 
nents of the droplet velocity, and the initial mass flow 
rate associated with each different droplet group as 


NASA/CR— 2012-217294 


20 



Table 6. Input file: nccJiquid_injector.in.01. 

Input file content 

comments 

heading 

title of a description of this file 

heading 

title of property 

nolc(nJ) 

denotes the total number of initial liquid components 
in the n i-tli injector. 

heading 

title of property 

out ^string 

chemical symbols of initial liquid species: e.g., C6H12. 

(ymki(n_i,n_l), 

nd=l,nolc(ni)) 

initial mass fractions of liquid species 

heading 

title of controlling parameters 

tdrop(n_i), 
isuperhd(n _i), 
ijetapr(nJ), 
vmfdrop(n_i) 

tdrop(n_i) denotes the initial droplet temperature. 

If isuperhd(nJ) = 1, it invokes the superheat vaporization 
model of Zuo et al. [62], and Schmehl et al. [63-64]; If isuperhd(n_i) 

= 2, it invokes the superheat vaporization model of Lee et al. [74]; & 
if isuperhd(n_i) = 0, it selects a non-superheat vaporization model. 

if ijetapr(mi) = .true., all the JP8 fluid properties are calculated by 
the computer coding provided by UTRC as described in Lee et al. [74]. 

vmfdrop(n_i) denotes the initial quality ( = gas mass fraction) of a 
superheated fluid. It only applies when ijetapr(nJ) = .true. 

heading 

title of controlling parameters 

atomization(nJ), 
drop_breakup_ 
model(nJ), 
spray_table(n_i), 
steady _spray_model 

If atomization(nj) = .true., the initial droplet conditions are 
determined based on the use of a primary atomization model. 
Otherwise they are determined from either known 
experimental conditions or a widely-used correlation. 

If drop-breakup_model(nJ) = .true., it invokes the secondary 
droplet breakup option. 

If spray .table = .true., initial droplet location, velocity, size, 
and mass flow rate are input through the file - ncc .liquid -table. in. 1. 
Otherwise they are determined by some other means. 

If steady .spray .model = .true., it invokes a steady state 
spray model commonly used in many spray codes, where 
whenever the spray solver is called, after first introducing 
a new group of spray particles, it continues with the 
liquid-phase computation until after all the particles are taken 
out of the computational domain (note: It is NOT 
recommended to use this option as a steady state calculation 
could be better arrived at by making use of some features of 
the unsteady option). The solution from the unsteady option 
(steady .spray unodel = .false.) is determined based on the 
values assigned to the controlling time steps - dtml, dtil, 

& dtgl, which are internal to the spray solver. 
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Table 6. Input file: nccJiquid_injector.in.01 (continued). 

heading 

title of controlling parameters 

no_of_holes (n_i) , 
no_of_streams(n_i), & 
no _of _drople t -groups (n J ) 

no_of_holes(nJ) = The number of holes per injector (note: 
when atomizationQ = .true, or spray .table = .true., the 
no of holes(n i) is set equal to 1). 


no _of .streams (n J) = The number of droplet streams per hole. This 
variable is introduced mainly to distinguish the variation in the 
droplet groups based on angular orientation (note: when atomizationQ 
.true, or spray .table = .true., the no_of_streams is set equal to 1). 

no_of_clroplet_groups(nJ) = The number of droplet groups per 
stream. This variable is introduced mainly to distinguish the 
droplet groups based on the droplet-size variation (note: 
when spray .table = .true, or atomizationQ = .true., all different 
injected droplets are lumped together into a single variable, 
no_of_droplet_groups() ) . 

heading 

title of controlling parameters 

lmdis(n_i), 

smdm(nJ), 

cone(nJ), 

size_min(nJ), 

size_max(n_i), & 

stochastic(n_i) 

If lmdis(nJ) = 1, it invokes the droplet size distribution as given 
by Eq. (72); If lmdis(nJ) = 2, it assumes a uniform droplet distribution 
between the maximum and minimum limits as specified by size_max() 

& size_min(); & If lmdis(nJ) = 3, it invokes the droplet size 
distribution as given by Eq. (73). 

smdm(nJ) = Sauter mean diameter 

If cone(nJ) = .true., it activates a 3D solid or hollow cone spray 
configuration as shown in Fig. 2. Otherwise it activates a 2D configuration 
depending on the value chosen for the logical variable - axisymmetric. If 
it is set equal to .true., it invokes the axis-of-symmetry case as shown 
in Fig. 3 otherwise it invokes the planar case as in Fig. 4. 

size_min(n_i) & size_max(n_i) = The variables are associated with 
the lmdis(nJ) = 2 option. 

stochastic (nJ) is used in conjunction with Option (2) of Sec. 9. If set 
equal to .true., it introduces some randomness into the determination 
of initial droplet velocity distribution. 

heading 

title of controlling parameters 

((x_inj(n_i,nx), 

yinj(n_i,nx), 

z_inj(nJ,nx), 

flowf(n_i,nx), 

v Jnj (n J,nx) , 

alpha_inj (n_i,nx) , 

beta_inj(nJ,nx), 

(x_inj(nJ,nx),y_inj(nJ,nx),zJnj(n_i,nx)) = spatial coordinates 
of the initial injector location. 

flowf(nJ,nx) = injector mass flow rate, (units - kgm/s for 3D & 
axis-of-symmetry & kgm/s/m for 2D planar. REMEMBER FOR AXIS-OF- 
SYMMETRY, IT IS THE TOTAL FLOW RATE OVER 360 DEGS.) 
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Table 6. Input file: nccJiquid_injector.in.01 (continued). 

theta jnj (nj,nx), 
dtheta jnj (n j,nx) , 
swlr_angle(n j,nx) ) , 
nx=l,no_oOioles(n j)) 

v_inj(nJ,nx) = droplet injection velocity, m/s. 
alpha Jnj (n_i,nx) = angle of rotation from the x-y plane, 
beta Jnj (nj,nx) = angle of rotation from the x-axis. 
theta jnj (nj,nx) = cone angle. 

dtheta jnj (n j,nx) = half-cone angle (note: Although dtheta jnj (n j,nx) 
= theta jnj (n j,nx)/2 for SOLID CONE SPRAY, it is invoked by 
setting dtheta jnj (nj,nx) either equal to 0 or theta jnj (nj,nx)/2) . Refer 
to Figs. 2-4, for a better understanding of the angular representation. 

swlr_angle(n j,nx) = swirl angle. 

heading 

title of controlling parameters 

(atom_type(n j,nx), 
breakup _type(n_i,nx), 
dia holefn i,nx), 
delp jnj (n_i,nx) , 
liq_vel(n j,nx), 
gas_u(n j,nx), 
gas_v(nJ,nx), 
gas_w(n j,nx), 
pcLstart(nJ,nx), 
pcl_end(nJ,nx)), 
nx= 1 , no _of Jroles (n j ) ) 

atom jype(n j,nx) = 0 => no primary atomization specified, 

= 1 => blob jet primary atomization model, 

= 2 => sheet breakup primary atomization model, 

= 3 => air blast primary atomization model, 

= 4 => BLS primary atomization model. 

breakup jype(n j,nx) = 0 => no secondary droplet breakup model specified, 
= 1 => Rayleigh-Taylor secondary droplet breakup model, 

= 2 => TAB secondary droplet breakup model, 

= 3 => ETAB secondary droplet breakup model. 

diaJrole(n j,nx) represents the injector orifice diameter (m). 

delp jnj (nj,nx) represents the pressure drop across the injector (Pa) 

(note: needed for pressure swirl atomization only). 

liq_vel(n j,nx) represents the liquid velocity from annulus (m/s) 

(note: needed for air-blast atomization only). 

gas_u(n j,nx), gas_v(n j,nx), gas_w(n j,nx): gas velocity components 
in annulus adjacent to liquid film (m/s) (note: needed for air-blast 
atomization only). 

pcl_start(n j,nx): starting parcel number for the annular air-blast injector. 
pcl_end(n j,nx)): ending parcel number for the annular air-blast injector. 
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Fig 2a. Geometrical details of fuel injection for 
a 3D solid or hollow cone spray. 



§ = azimutal angle, deg 
number of rays = 8 ((|)=45) 
For a 3D cone, 
no_o f_st reams ( ) should 
defined as a multiple of 8 
number along each ray = 
no_of_st reams 0/8. 


Fig 2b. Initial spray particle orientation in a 
circular cross-section. 
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9 = cone angle, deg 
d0= half-cone angle 



Fig 3a. Geometrical details of fuel injection for 
an axis-of-symmetry case. 


• 00 

* 

• 01 


For axis of symmetry: 
no_of_streams ( ) is distributed 
between 01 to 00 . 


; C 


Fig 3b. Initial spray particle concentration in an 
axis-of-symmetry case. 
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9 = cone angle, deg 



Fig 4a. Geometrical details of fuel injection for 
a 2D planar case. 


• 00 

♦ 01 For 2D : 


no. 

_of_streams ( ) 

i / 2 

is between 

01 

to 00, and 



no. 

_of_streams ( ) 

/ 2 

is between 

LO 

to LI. 




* LI 

V 

* 

'• LO 


Fig 4b. Initial spray particle concentration in a 
2D planar case. 
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Table 7. Input file: nccj3pray_table.in.01. 


Input file content 

comments 

nos(n_i) 

denotes the total number of droplet groups in the n j-th 
injector. 

( ni, xx Jnj, yyJnj, 

ni = number of the droplet group. Its value ranges between 

zz _inj ,uu_inj , vv_inj , 

1 to nos(nj). 

ww jnj ,r jnj , flcLd) 
,ni=l,nos(n j)) 

(xx inj.yy inj,zz inj) = spatial coordinates. 


(uu jnj,vv jnj,ww jnj) = velocity components. 


r jnj = droplet size in radius. 


flcLd = mass flow rate of the ni-th droplet group (note: 
SUMMATION OF fld_d OVER ALL THE nos(n j) 
DROPLET GROUPS IS EQUAL TO THE TOTAL MASS 
FLOW RATE OF THE INJECTOR). 


described in ncc_spray .table. 1 of Table 7. This ta- 
ble provides a complete description of the following 
variables: nos(nj), (ni,xx_inj, yyJnj, zzdnj, uujnj, 
vvJnj, wwTnj, rJnj, flcLd), ,ni=l,nos(n_i)). The ini- 
tial inputs specified through a spray table should be 
representative of the integrated averages of the ex- 
perimental conditions [10-13]. 

Option (2): 

In this option, we need to specify several pa- 
rameters including no_ofJroles(), no_of_streams(), & 
no_of_droplet_groups(). And depending on what is 
specified for the input of the integer variable, lmdis, in 
ncc jnjector.in.l, the droplet-size distribution within 
each one of the streams is determined based on the 
following three choices: 

• In one choice, it is calculated based on a corre- 
lation typical of those widely used in describing 
the initial droplet size distribution [4] : 

- 16 - 98 ( 4 ) 04 ^ 

^32 

(72) 

where n is the total number of droplets and dn is 
the number of droplets in the size range between 
d and d + dd. This correlation also requires the 
specification of Sauter mean diameter, ^32. Fig. 
5 shows the droplet size distribution generated 



Droplet dldmiir, n I crons 


Figure 5 Droplet-size distribution. 

by Eq. (72) for a case studied in [10]. The solid 
line shows the droplet number variation versus 
drop size and the dashed line shows the inte- 
grated mass variation with drop size. The drop 
size distribution within the spray is represented 
by a finite number of droplet classes as given by 
the variable, no_of_droplet_groups(). 

• In the second choice, the initial droplet sizes are 
distributed evenly between two specified maxi- 
mum and minimum drop sizes - size_max() & 
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size_min() as specified in Table 6. And the size 
interval depends on the value specified for the 
variable, no_of_droplet_groups(). However, this 
option is applicable only in certain special cases. 

• In the third choice, the initial droplet size dis- 
tribution is calculated by making use of a cliped 
probability distribution function as given by 

P(x) = [ f{x)dx/ [ f(x)dx (73) 
Jo Jo 

where /(x) = (1.0 — exp(— x){l+x+\x 2 + \x^). 
After a random sampling of 0 < % < 12.0, the 
initial droplet diameter is determined as follows: 
(0.3 < 4P(y) < 1.5)d32- If the value for 4P(x) 
goes beyond the limits of 0.3 and 1.5, the ran- 
dom sampling is reiterated until that number 
falls within the range of the lower and upper 
bounds. 

Depending on what is specified for the logical 
variable, cone, of Table 6, the droplet velocity dis- 
tribution amongst various streams of a given hole is 
calculated by assuming the spray to be either a solid 
or hollow cone spray. A graphical illustration of three 
different cone configurations are shown in Figs. 2 to 
4. Figs. 2a & 2b refer to a 3D case, Figs. 3a & 3b 
to an axisymmetric case, and Figs. 4a & 4b to a 2D 
planar case. It is noteworthy that in an axisymmet- 
ric case, the x-axis is assumed to be aligned with the 
axis-of-symmetry. 

Fig. 2a shows the geometric details of a hollow 
cone spray in 3D where 9 is known as the cone angle, 
dO is the half-cone angle, and for a solid cone spray 
d9 = 9/2. And a represents the angular rotation 
from the x-y plane and (3 is the angular rotation from 
the x-axis. Fig. 2b shows initial spray stream orien- 
tation in a circular cross section. We try to simulate 
the spray by a finite number of streams as given by 
the variable, no_of_streams(). Each one of the dark 
circles in Fig. 2b represents a different stream. The 
streams are distributed evenly along each one of the 
different rays which are separated from each other 
with an angle of separation as given by 4>. Currently, 
we have hard-coded the angle of separation to be 45° 
which means we have restricted the number of rays 
to be eight in 3D. Therefore, when specifying a num- 
ber for no_of_streams(), it should be borne in mind 
that it should be a multiple of eight. And the num- 
ber of streams along each one of the rays, therefore, 
becomes equal to no_of_streams()/8. In Fig. 2b, the 


no_of_streams() has a value of 32 and, therefore, the 
number of streams along each one of the rays for this 
case is 4. 

Fig. 3a. shows the geometric details for an 
axisymmetric case. Since the computations are per- 
formed only in the first quadrant of the x-y plane, all 
of the specified number of streams, no_of_streams(), 
are distributed evenly between the lower and upper 
halfs as shown in Fig. 3a and 3b. Here, the upper half 
refers to the region between 01 to OO and the lower 
half refers to LO to LI. It is noteworthy that each 
droplet group in a given stream represents a circular 
ring of liquid for an axisymmetric case. 

The geometric details for a 2D case are shown 
in Figs. 4a & 4b. Here, it is noteworthy that 
each droplet group in a given stream represents a 
planar sheet of liquid. All the specified number of 
streams, no_of_streams(), are distributed evenly over 
both sides of the cone center as shown in Fig. 4b. 

For each one of the different injector holes, 
no_ofJroles(), it requires the specification of the fol- 
lowing parameters as described in ncc -injector. in. 1 
of Table 6 (However, because of the geometric dif- 
ferences that exist between 3D, 2D, and axisymmet- 
ric configurations, some of the input parameters may 
have different units as noted below): 

• The initial (x, y, z) coordinates of the hole loca- 
tion. 

• The mass flow rate per hole - however, the defi- 
nition of the units for the injector mass flow rate 
per hole differs: it is kgm/s for 3D & the axis- 
of-symmetry and kgm/s/m for 2D planar (note: 
In an axis-of-symmetry, the specified mass flow 
rate per hole refers to the entire mass flow rate 
over 360 degrees). 

• The following variables define the angular ori- 
entation: cti n j =angle of rotation from the x-y 
plane, /3i n j = angle of rotation from the x-axis, 
9 in j = cone angle, & dO in j = half-cone angle 
(note: Although d9i n j = 9i n j / 2 for a solid cone 
spray, a specified value of zero for d9i n j also in- 
vokes a solid-cone spray configuration). 

• The variable, swlr_angle(), allows a means to 
specify the tangential component in the case of 
both 3D and axisymmetric sprays. 

Stochastic Option 
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In the above description of Option (2), the ini- 
tial droplet conditions are determined based on a de- 
terministic approach. However, such an approach 
when employed in 3D may lead to time consuming 
calculations requiring the use of many droplet groups. 
With this in mind, we wanted to introduce some ran- 
domness into the calculation of the initial droplet ve- 
locity distribution as an option. This option can be 
invoked by setting the logical variable, stochastic(), 
to be .true.. In a stochastic approach, the initial an- 
gular orientation of a particle could be randomized in 
a number of different ways. The obvious way would 
be to determine the initial droplet direction based on 
a complete randomization in a field as defined by the 
geometric parameters: solid cone angle, 9 , half cone 
angle, d6 and/or azimuthal angle, </>. Another way 
is try to retain the basic features of the determinis- 
tic approach where the initial droplet directions are 
taken as the means in a random field with the ran- 
domness localized to within the surrounding angular 
sub-segment of the respective droplet as represented 
by 86, Sdd , 868<j), or 8d68<j). We have chosen to go 
with the later option as it preserves the same basic 
structure but provides some benefits of the stochastic 
approach. 

Option (3): 

Here many aspects of fuel injection remain the 
same as in Option (2) but the initial droplet size 
and velocity distribution may differ depending upon 
the primary atomization model selected. More de- 
tails on the available primary atomization models 
can be found in the appendix. For each one of the 
different injector holes, no.ofJrolesQ it requires the 
specification of all the parameters as in Option (2). 
However, it doesn’t make use of the values set for 
size_max() & size_min(). Also, it requires the spec- 
ification of some additional parameters as described 
in ncc_injector.in.l of Table 6. However, there exist 
some differences with regard to their usage depending 
upon the primary atomization model. Some major 
differences are highlighted below: 

• The integer variable, atom_type() of Table 6, is 
used to select the primary atomization model of 
one’s choice. 

• dia_hole() needs to be specified for all primary 
atomization models. 

• delp_inj() needs to be specified only with the 
sheet breakup primary atomization model and 
it is used in Eq. (9) of the appendix. 


• liq_vel(nJ,nx), gas_u(n_i,nx), 

gas_v(nJ,nx), gas_w(nJ,nx), pcl_start(n j,nx), & 
pcl_end(nJ,nx)) are required only for the air- 
blast primary atomization model. 

• In all of the primary atomization models except 
for the air-blast, the total number of droplet 
groups (particles) to be injected at the time 
of every injection period is specified though 
the use of the variable, no_of_droplet_groups(). 
For the air-blast, it is given by pci end(n i,nxj- 
pcl_start (n J,nx) + 1 . 

• In this option, all of the injected droplet groups 
as specified by either no_of_droplet_groups() 
or pcLend(nJ,nx)-pcLstart(nJ,nx)+l are dis- 
tributed randomly across the entire cross-section 
of any selected 2D, axisymmetric, or 3D spray 
configuration. So it is recommended that 
no_of_streams() be set equal to one as no dis- 
tinction is made between streams. 

The secondary droplet breakup option can be 
invoked regardless of how the initial droplet condi- 
tions are generated using Options (1) to (3). The sec- 
ondary droplet breakup option can be invoked by set- 
ting the logical variable, drop .breakup _model(), to be 
.true.. And then the integer variable, breakup_type() 
of Table 6, is used to select the secondary droplet 
breakup model of one’s choice. More details on the 
available secondary droplet breakup models can be 
found in the appendix. 

10 SPRAY SOLUTION ALGORITHM 

In order to evaluate the initial conditions that 
are needed in the integration of the liquid-phase equa- 
tions, we first need to know the surrounding gas- 
phase properties at each particle location. But in or- 
der to evaluate the gas-phase properties, it is first nec- 
essary to identify the computational cell in which a 
given particle is located. It is a trivial task to track a 
particle in the regular rectangular coordinates. How- 
ever, the particle tracking becomes complicated when 
the computational cell is no longer rectangular in the 
physical domain and it becomes even more compli- 
cated when the particle search is performed within 
the context of parallel computing. 

We have developed and implemented an effi- 
cient particle-tracking algorithm for use with parallel 
computing in an unstructured grid. The search is ini- 
tiated in the form of a local search originating from 
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Fig. 6 A vector illustration used in the particle search 
analysis. 

the computational cell in which the same particle was 
found to be located in the previous time step. The 
location of the computational cell is then determined 
by first evaluating the dot product of x pc . a n = \x pc \ 
\a n \ cos ((/)), where x pc is the vector defined by the 
distance between the particle location and the center 
of the n-face of the computational cell, a n is the out- 
ward area normal of the n-face as shown in Fig. 6, 
and (f) is the angle between the two vectors. 

A simple test for the particle location requires 
that the dot product be negative over each and ev- 
ery one of the n-faces of the computational cell. If 
the test fails, the particle search is carried over to the 
adjacent cells of those faces for which the dot prod- 
uct turns out to be positive. Some of those n-faces 
might represent the boundaries of the computational 
domain while the others represent the interfaces be- 
tween two adjoining interior cells. The search is first 
carried over to the adjacent interior cells in the di- 
rection pointed out by the positive sign of the dot 
products. The boundary conditions are only imple- 
mented after making sure that all other remaining 
possibilities point towards a search exterior of the 
computational domain. This implementation ensures 
against any inadvertent application of the boundary 
conditions before the correct interior cell could be 
identified. 

After the gas-phase properties at the particle 
location are known, the solution for the ordinary dif- 
ferential equations of particle position, size, and ve- 
locity are advanced by making use of a second-order 
accurate Runge-Kutta method. The partial differen- 
tial equations governing the droplet internal thermal 
and mass transport are integrated by making use of 
a fully implicit Newton-Raphson iteration method. 

Finally, the liquid-phase source terms of the gas- 
phase conservation equations (1-4) are evaluated by 


making use of a time-averaging method. 

11 DESCRIPTION OF COMPUTER 
FLOWCHART, & THE TIME- AVERAGING 
SCHEME USED IN THE GAS-PHASE 
SOURCE-TERM CALCULATIONS 

• In order to know more about the time-averaging 
method, we need to know first about the three 
different time steps that are internal to the spray 
code: A t m i, A tu, and Af g /. 

A t m i - the actual time step used in integrating 
the liquid-phase equations which is determined 
based on the smallest of the different time scales 
associated with various rate-controlling phenom- 
ena of a rapidly vaporizing droplet. 

A tu - the injection time step. It is the time 
step at which a new discretized parcel of different 
droplet groups are introduced into the computa- 
tion. 

A t g i - the global time step. Its introduction 
seems to provide better convergence in both un- 
steady and steady-state computations. 

• When the spray solver is called it advances the 
liquid phase equations over a number of itera- 
tions as determined by the ratio of Atgi/At m i. 

• It then evaluates the time-averaged contribution 
of the liquid-phase source terms, S g i, of the gas- 
phase governing equations (1-4) as follows: 

M Ay- 

®- - £ < 7 “> 

m — 1 y 

where 


M 

^rnl = A tgl (75) 

m — 1 

• The values for A t m i, A t g i, & A tu are specified 
in the input file, nccdiquid_solver.in, of Table 8. 

• In steady-state computations, it is recommended 
to use for both At g i and A tu a value of about 
1 ms which is roughly equivalent to the aver- 
age lifetime of the droplets for a typical reacting 
spray encountered in conventional low-pressure 
gas-turbine combustors. 
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Table 8. Input file: nccJiquid_solver.in. 

Input file content 

comments 

heading 

title of controlling parameters 

ldread,ispray_mod, 

( when_st art _spr ay (n j ) , 
n j=l,no_of_injectors) 

If ldread = .true., restarts the calculation from the data stored from 
the previous computation. Otherwise initiates a new spray computation. 

ispray_mod controls the calls to the spray solver. The spray solver is called 
once at every other CFD iteration as given by the number, ispraycmod. 

when_start_spray() represents the CFD iteration where the computations for 
the n_ith injector are initiated. 

heading 

title of controlling parameters 

dtml, dtgl, & 
dtil 

dtml = time step for advancing the liquid phase equations. 

dtgl = global time step. Whenever the spray solver is called, it advances 
the spray computations over a period of dtgl before returning control 
over to the calling routine. To be more precise, it advances the 
liquid phase equations over a number of time steps as determined by 
(dtgl/dtml). 

dtil = injection time step. It is where a new group of droplets are 
introduced into the computation. 


The averaging scheme could be explained better 
through the use of a flow chart shown in Fig. 7. 
The main spray solver is invoked with a controlling 
routine, DCLR, which, then, executes the following 
steps: 

1. It first initializes the source terms to zero. 

2. Updates the global time, t g i, based on A t g i- 

3. Checks to see if t m i < t g i < t m i + A t m i. If it 
is, it returns control over to the calling routine 
and supplies the other flow solvers, e.g., flow or 
EUPDF, with the source terms, S g i, of Eqs. (1)- 
(4). If not, it proceeds with the next step. 

4. Checks to see if it it is time to introduce a new 
group of particles. 

5. Proceeds with solving the liquid-phase equations 
with calls to the following routines: 

• Calls the particle tracking routine and as- 
signs particles based on the parallel strategy 
implemented. 

• Interpolates gas-phase properties at the 
particle location. 


• Advances liquid-phase equations and, then, 
deletes any particles that are no longer 
needed in the computations. 

6. Evaluates the liquid-phase source-term contribu- 
tions, S rn i , of Eq. 74. 

7. Updates the time, t m i, based on At m i. 

8. It then goes back to step (3) and repeats 
the whole process again until the computa- 
tions are completed over a global time step: 
£m=i At ™l = A V 


12 IMPLEMENTATION OF BOUNDARY 
CONDITIONS 

The spray code supports most of the boundary 
conditions that are in use in the current version of 
the NCC CFD module but not all. Here, we would 
like to highlight some of the boundary conditions: 

• In case of the droplet impingement with the com- 
bustor walls, the droplets may evaporate, move 
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Fig. 7 The flowchart of LSPRAY-IV. 
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along the wall surfaces, and/or reflect with re- 
duced momentum. The physics of droplet im- 
pingement with the walls is not completely un- 
derstood. In our present calculations, it is as- 
sumed that the droplets, after having lost most 
of their momentum upon impingement with the 
walls, move along the wall surfaces with a veloc- 
ity equal to that of the surrounding gas. 

• The implementation of the periodic boundary 
conditions becomes rather complicated as it re- 
quires taking into consideration several aspects 
arising from a particle leaving the computational 
domain from one periodic boundary has to reen- 
ter the domain back from a corresponding sec- 
ond periodic boundary with appropriate condi- 
tions. The periodic boundary conditions are im- 
plemented with the help of some appropriately 
defined transformation matrices. 

• The symmetric boundary condition is imple- 
mented in such a way to satisfy the criterion that 
for every particle crosssing the symmetry line, a 
similar one re-enters the domain in a direction 
given by the reflection off of the symmetry line. 

• When the particles try to move out of the exit 
boundary, they are removed out of the compu- 
tation. 


13 DETAILS OF COUPLING BETWEEN 
LSPRAY-IV AND OTHER SOLVERS 

• The spray code is designed to be a stand-alone 
module such that it could easily be coupled with 
any other unstructured-grid CFD solver and the 
same holds true for EUPDF. (However, some 
grid-related parameters such as area vectors, grid 
connectivity parameters, etc. need to be sup- 
plied separately). 

• The spray solver needs information on the gas 
velocity and scalar fields from the other solvers 
and, then, it in turn supplies the liquid-phase 
source terms. 

• The PDF solver needs information on the mean 
gas velocity, turbulent diffusivity and frequency 
from the CFD solver and the liquid-phase source 
terms from the spray solver, and then it in turn 
provides the solution for the scalar (species and 
energy) fields to the flow and spray solvers. 


• It should also be noted that both the PDF and 
spray solvers are called once at every other spec- 
ified number of CFD iterations. 

• All of the three solvers (LSPRAY-IV, EUPDF, 
and CFD) are advanced sequentially in an iter- 
ative manner until a converged solution is ob- 
tained. 

• All three codes (EUPDF, CFD, and LSPRAY- 
IV) were coupled and parallelized in such a way 
to achieve maximum efficiency. 

The coupling issues could be better understood 
through the use of a flow chart shown in Fig. 8. 
It shows the overall flow structure of the combined 
CFD, LSPRAY-IV, and EUPDF modules. Both the 
PDF and spray codes are loosely coupled with the 
CFD code. The spray code is designed in such a 
way that only a minimal amount of effort is needed 
for its coupling with the flow and PDF solvers. The 
present version of the spray module relies entirely on 
the use of Fortran common blocks for its informa- 
tion exchange with other modules. Even this reliance 
should entail only few changes to be made within the 
spray code for its linkage with different solvers. The 
PDF code is also structured along similar coupling 
principles. 

The flow chart of Fig. 8 contains several blocks 
- some shown in black and/or solid lines and the oth- 
ers in color and/or dashed lines. The ones in solid 
blocks represent the flow chart that is typical of a 
CFD solver. The ones in dashed blocks represent the 
additions arising from the coupling of the spray and 
PDF solvers. 

The coupling starts with the calling of the 
subroutine - spray pdf read parameters, which 
then reads the spray control parameters from the 
input file, nccJiquidsolver.in of Table 8. This 
table provides a detailed description of the fol- 
lowing input file variables - ldread, ispray_mod, 
(when_start_spray(n _i), nJ=l,no_ofJnjectors), dtml, 
dtgl, & dtil. The coupling is then followed by the 
calling of the pdf_int_rerun subroutine. It ini- 
tializes PDF computations and also it may restart 
the PDF computations if needed from the data 
stored from a previous iteration. Similarly, we 
call spray int rerun for the spray computations. 
It is noteworthy that the spray computations can 
be restarted by reading the data from the restart 
files: nccJiquidjpai'ams.out , nccJiquidjresults.db, 

& nccJiquid_resultsjmc.db. The restart capability 
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SPAWN MULTIPLE PEs FOR PARALLEL COMPUTING 

i 

READ PARAMETERS & GEOMETRIC DATA 

, -L- 

I CALL SPRAY PDF READ PARAMETERS I 


INITIALIZE 
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READ PDF RESTART FILES 
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Fig. 8 The overall flowchart of the combined CFD, LSPRAY-IV, and EUPDF solvers. 
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is invoked by setting the logical variable, ldread, of 
the input file, nccJtiquidsolver.in of Table 8, to be 
true. Otherwise, the spray computations are initial- 
ized to start from the beginning. Then, the coupling 
proceeds with the calling of the following subroutines: 
dclr for integrating the spray calculations and eupdf 
for the Monte Carlo PDF. The input variable, is- 
pray_mocl of Table 8, controls the calls to the spray 
integrating routine. The spray solver is called once at 
every other number of CFD iterations as specified by 
ispray_mocl. And the first call to the spray solver is 
controlled by the input variable, when_start_spray() 
of Table 8, which represents the starting CFD itera- 
tion number from which the spray computations are 
initiated. Finally, the coupling ends with the call- 
ing of a subroutine, spray pdf output, which will 
create a set of new restart files. 

14 DETAILS OF FORTRAN 
SUBROUTINES & FUNCTIONS 

Table 9 provides a list of all the Fortran sub- 
routines developed as a part of the spray module. 
This table also provides information on all the For- 
tran functions. It also describes the purpose of all 
the individual subroutines and functions. Similary, 
Table 10 provides a list of all the Fortran functions 
and subroutines associated with the implementation 
of the high pressure equation of state as described in 
Sec. 3, and the high pressure behavior on the gas- 
phase transport properties as in Sec. 4. 

15 DETAILS OF PARALLELIZATION 

There are several issues associated with the par- 
allelization of both the spray & PDF computations. 
The goal of the parallel implementation is to extract 
maximum parallelism so as to minimize the execu- 
tion time for a given application on a specified num- 
ber of processors [37]. Several types of overhead costs 
are associated with parallel implementation which in- 
clude data dependency, communication, load imbal- 
ance, arithmetic, and memory overheads. The term 
arithmetic overhead is the extra arithmetic opera- 
tions required by the parallel implementation. Mem- 
ory overhead refers to the extra memory needed. Ex- 
cessive memory overhead reduces the size of a prob- 
lem that can be run on a given system and the 
other overheads result in performance degradation 
[37]. Any given application usually consists of sev- 
eral different phases that must be performed in cer- 
tain sequential order. The degree of parallelism and 


data dependencies associated with each of the sub- 
tasks can vary widely [37]. The goal is to achieve 
maximum efficiency with a reasonable programming 
effort [37]. 

In our earlier work, we discussed the parallel 
implementation of a spray algorithm developed for 
the structured grid calculations on a Cray T3D [10]. 
These computations were performed in conjunction 
the Monte Carlo PDF method. The parallel algo- 
rithm made use of the shared memory constructs ex- 
clusive to Cray MPP (Massively Parallel Processing) 
Fortran and the computations showed a reasonable 
degree of parallel performance when they were per- 
formed on a NASA LeRC Cray T3D with the number 
of processors ranging between 8 to 32 [10]. Later on, 
the extension of this method to unstructured grids 
and parallel computing written in Fortran 77 with 
PVM or MPI calls was reported in [11-15]. The lat- 
est version in Fortran 77/90 offers greater computer 
platform independence. In this section, we only high- 
light some important aspects of parallelization from 
Refs. [10-15], 

• Both the EUPDF and CFD modules are well 
suited for parallel implementation. For the gas- 
phase computations, the domain of computation 
is simply divided into n-Parts of nearly equal 
size and each part is solved by a different pro- 
cessor. Fig. 9a illustrates a simple example of 
the domain decomposition strategy adopted for 
the gas-phase computations where the total do- 
main is simply divided equally amongst the avail- 
able computer processing elements (PEs). In 
this case, we assumed the number of available 
PEs to be equal to four. 

• But the spray computations are more difficult to 
parallelize for the reasons summarized below: 

(1) Non-uniform nature of spray distribution: 
Most of the particles are usually confined to a 
small region where the atomizer is located. 

(2) Dynamic nature of Lagrangian particles: 
Particles keep moving between different sub- 
divided domains of an Eulerian grid (grid used 
in the CFD computation) which are assigned to 
different processors. While some new particles 
are introduced at the time of fuel injection, some 
others are taken out of computation. 

Conceptually, there are several ways to paral- 
lelize the spray computations, we, however, devel- 
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Table 9. Description of LSPRAY-IV Fortran subroutines & functions. 

Function 

Purpose of the Function 

blasiu(x) 

This function returns a solution for the function, f(Bk), 
of Eq. (21) for use in computing the droplet regression rate. 

Subroutine 

Purpose of the Subroutines 

axisym_part _adj ust 

This routine was written by Jeff Moder and essentially 
performs the same functions as the intpla routine but 
modified for application with a 2d-axisymmetric case. 

calc_F_suvw() 

It evaluates the right-hand sides of ODEs for both 
droplet size and velocities, and then provides the 
solution to the calling routine, chaslv. 

chaslv 

This routine has two main functions: 

(1) It integrates the liquid phase equations. 

(2) It removes the particles that are no longer needed 
in the computation. 

dclr 

This routine is called once at every other CFD iteration as 
specified by ispray_mod. It is primarily a controlling 
routine for spray computations. 

This routine has the following functions: 

(1) It initializes the source terms to zero. 

(2) Checks to see if new particles need to be introduced. 

(3) Advances liquid phase equations over an allowable or 
pre-specified time step, dtml, with calls to the following 
routines: 

intplal - Interpolates the gas phase properties at the particle 
location. 

chaslv - Advances liquid phase equations. 

intpla - Identifies computational cells and PEs associated 
with spray particles. 

sprips - Evaluates the liquid phase source term contributions 
of the CFD and PDF equations. 

(4) (a) For unsteady spray computations (steady_spray_model=F), 
this routine is called once at the beginning of every 

global time-step, dtgl. This is accomplished by continuing 
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Table 9. Description of LSPRAY-IV Fortran subroutines & functions (continued). 

Subroutine 

Purpose of the Subroutines 


with steps (2) and (3) until the computations are completed 
over one single global time step, cltgl. (b) For steady spray 
spray computations (steady .spray _model=T), this routine is 
only called once to solve entire lifetime of all particles 
injected. It repeats step (3) until there are no more spray 
particles left to be integrated (note: The particles are taken 
of computation when they reach either a certain size of 
negligible proportion or exit out of the computational domain). 

(5) Returns control over to other solvers, e.g. CFD and EUPDF, 
and supply them with the source terms averaged over dtgl. 

dropdis(rhol, 

flowdum,sr, 

fld,smd,nofg) 

This routine provides initial droplet distribution from 
the following correlation (also, appears as Eq. 72): 

dn/n = a((D/D 32 ) al P)exp(-b((D/D 32 ) bet ))dD/D 32 
where a, b, alp, and bet are constants. 

drop jc(n_i, 
nmih,nmis, 
nmip) 

It is called by the spray jnainJnjection routine in conjunction 
with spray _table(n_i))=F at the beginning of every injection 
time step, dtil, for introducing a new group of spray particles. 

It applies for the introduction of new particles in conjunction with 
Options (2) or (3) of Sec. 9. In Opt. (2) the initial conditions 
are specified from the known correlations and in Opt. (3) from 
the available primary atomization models. 

find_transport _ds ( 
ijle,chem_tnodel, 
number _of_species,hp_eos, 
temp gas, pres gas.y gas, 
rhogm,visgm,congm,difgm) 

It computes the following properties of a gas mixture 
at the droplet interface by making use of the one-third 
rule of Eq. (46): molecular viscosity, gas density, 
thermal conductivity, and diffusion coefficient. 

find_xyzface 

(i,xcfac,ycfac, 

zcfac) 

This routine computes x, y, and z locations of all the face 
centers of an element, i. This information is used in the 
particle search algorithm. 

get _liq_t at .properties 
(tmp,pmp,j) 

It computes the following variable properties of a liquid 
mixture: density, specific heat, molecular viscosity, 
gas density, thermal conductivity, and diffusion coefficient. 
This routine also computes the surface tension of a 
multi-component liquid fuel. 

get _liq_t at .properties 1 

(tmp,pmp,boil_tem, 

cp_sph,hvap_sph,icuin) 

This subroutine is used in conjunction with the superheat 
vaporization model and computes the following variable 
properties of the liquid fuel: boiling temperature, specific 
heat at constant pressure, surface tension, heat of vaporiza- 
tion, molecular viscosity, & liquid density and volume. 

get_liq_tat_properties_sp 

(tmpl,tmp2,cpltl, 

cplt2,j) 

This subroutine is used in conjunction with the superheat 
vaporization model and computes the following variable 
properties of the liquid fuel: specific heat at droplet 
temperature, & specific heat at boiling temperature. 
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Table 9. Description of LSPRAY-IV Fortran subroutines & functions (continued). 

Subroutine 

Purpose of the Subroutines 

init_psi2_dsd 

This routine provides initial droplet distribution from 
a clipped probability distribution function as given by 
Eq. 73. 

intpla 

This routine performs three main functions: 

(1) Particle Tracking - It identifies the computational cell 
in which a particle is located. In parallel computing, 

it also means identifying the corresponding processor of the 
computational domain in which a particle is located. 

(2) It implements appropriate boundary conditions. 

(3) It reassigns the particles between different PEs based on 
the parallel strategy employed. 

intplal 

This routine interpolates the gas-phase properties at the 
particle location. In the present case, a simple first-order 
interpolation is employed. 

mimd_spray 

For computational elements where the neighboring cells are 
partitioned between different processors, this routine initializes 
the arrays, ipr_fr_id() and ile_fr_id(), by storing the information 
on the corresponding processor and element ID numbers of 
all the neighboring cells. This information is needed in order to 
track the particle movement between neighboring cells. 

mimd_spray_recv 

(i_recvfrom) 

This subroutine is called by mimd_spray in order to gather 
some relevant information from the neighboring processors. 

mimd_spray_sencl 

(i_sendto) 

This subroutine is called by mimd spray in order to send 
some relevant information to the neighboring processors. 

par_loc(xparz, 

yparz,zparz, 

ipare,iparp) 

Given the x,y,z coordinates of a particle location, the 
algorithm identifies the corresponding computational cell 
in which a particle is located. It also identifies the 
corresponding processor of the computational grid in which 
a particle is located. This information is used to locate 
newly injected droplet groups. 

prnspr 

It provides written output of the spray data in a standard format. 

spray _int_rerun 

This routine has two functions: (1) It initializes the spray 
solver based on the data read from various input files. 

(2) If it is a rerun, it restarts the computations from the 
stored data of a previous computation. 
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Table 9. Description of LSPRAY-IV Fortran subroutines & functions (continued). 

Subroutine 

Purpose of the Subroutines 

spray .main JnjectionQ 

This routine is called to inject a new group of spray 
particles for each injector at the beginning of every 
injection time step, dtil. 

spray .output 

This routine is called to create appropriate restart files 
for the spray computations and also to provide for debugging 
purposes written output of some useful gas and spray data in 
a standard format. 

spray _plot_output 

It provides written output of some useful spray data in 
a standard format. It can be used in conjunction with 
steady jspray .model = F. 

spray _read_parameters 
(ipid,ncells)) 

This routine is called once in the beginning to read and 
initialize some of controlling parameters associated with 
the spray computations. 

spray .read-transport 
( (ipid, fid, filename) 

This routine is called once to read liquid thermo and 
tranport data: various constants used in the evaluation 
the polynomials in temperature for constant pressure 
specific heat, thermal conductivity, and viscosity. 

sprips 

This routine is called to provide the liquid-phase source terms of 
Eqs. (l)-(4) for use in both CFD and Monte Carlo PDF solvers: 

smlc(i) = liquid-phase contribution of Eq. (1) of Section 
2. 

smlmx(i), smlmy(i), smlmz(i) = liquid-phase contribution 
of Eq. (3) of Section 2. 

smle(i) = liquid-phase contribution of Eq. (4) of Section 2. 

sy(il,iu,bb, 

clcl,aa,cc) 

It is a tri-diagonal matrix solver. It is used in the solution 
of both Eqs. (27) & (37). 

uvw_par(swlr_angle(nx) , 
angle _work, nx, ny, nn , 
v Jnj ,nmip,t_rotation, 
cone, n_cone_r ay s , 
cone .rotation, uloc, 
vloc,wloc,nJ) 

It computes the initial particle injection velocity for different 
geometric configurations of spray as shown in Figs. 2 to 4. 
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Table 10. Description of PENG_ROB Fortran subroutines & functions. 

Function 

Purpose of the Function 

p eng _rob .density 
(chem_model,ycomp, 
tij, pressure, 
number _of_species) 

It provides a solution for the density of a multi-component 
mixture based on the Peng-Robinson equation of state. 

peng_rob .pressure 
(chem_model,ycomp, 
tij , density, 
number _of_species) 

It provides a solution for the pressure of a multi-component 
mixture based on the Peng-Robinson equation of state. 

Subroutine 

Purpose of the Subroutines 

hp.diff.table 

(ct,pt,dabcvi) 

This subroutine provides a value for {Dab P)/(D, 4 b P) + 
based on the Takahashi high-pressure correlation for the 
calculation of the diffusion coefficients in a gas mixture. 

It is based on the interpolation from the data derived from 
a graphical illustration of Reid et al [16], “The Properties 
Gases and Liquids.” The retrieved data is stored in 
hp_diff_table.dat as a table of f(tpr_diff,ppr_diff) values. 

hp .diffusion 

(tij,pij,its, 

nwarn,diin_hp) 

It is called to evaluate the gas-phase mixture-averaged 
diffusion coefficients valid at high pressures. 

hp.transport 
(calc_name,its, 
tij,pij,nwarn, 
avisi Jrp , condi _hp ) 

It is called to evaluate the gas-phase transport properties 
valid at high pressures: viscosity and conductivity. 

p eng .rob .all _rhos ( 
chem_model,ycomp, 
tij, pressure, 
number _of .species , 
rhol, rho2, rho3, 
vvl, vv2, vv3, 
zl, z2, z3, 
aprgm,bprgm) 

It provides the roots for the cubic equation, Eq. (6), involving 
the compressibility factor, Z. One of the three roots of Eq. (6) 
yields Z v and the other Z ;. Jeff Moder made significant 
improvements to this routine. 

peng_rob_gen 

(xsp,ppr,tpr, 

its,rho0) 

It takes the following as the input variables: xsp() - the mass 
fraction of the species, ppr - pressure, tpr - temperature, and 
its - the number of species. And it returns as the output, rhoO - 
density based on the solution of the Peng-Robinson EOS. 

peng.rob_zcrit.rhor 

(ppr,tpr,its) 

It provides a solution for Zj C and p lr for ith species based 
on the solution of the Peng-Robinson EOS for a given value 
of P.; r and Tj r . It also requires the specification of the 
total number of species, its. 
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Fig. 9a An illustration of the parallelization 
strategy employed in the gas flow comp- 
utations . 
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Fig. 9b An illustration of the parallelization 
strategy employed in the spray computa- 
tions . 
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Table 11. CPU time (sec) per cycle 

versus number of PEs. 



Number of processors 

Solver 

Characteristic 

2 

5 

10 

CFD 

5 steps/cycle 

2.50 

1.25 

0.75 

EUPDF 

1 step/cycle 

6.5 

2.9 

1.9 

LSPRAY 

100 steps/cycle 

1.70 

0.64 

0.53 


Table 12. CPU Time (Sec) Per Cycle Versus Number of PEs. 


Number of Processors 

Solver 

Characteristic 

1 

2 

4 

8 

16 

Spray 

100 steps/cycle 

6.83 

5.29 

2.94 

1.64 

0.87 

Max. Spray Particles in a PE 

2695 

2097 

1165 

623 

312 

Min. Spray Particles in a PE 

2695 

598 

118 

14 

0 


oped and tested two different domain decomposition 
strategies [10-14]. 

• Strategy I: 

The Lagrangian particles were assigned fairly 
uniform amongst the available processors but the 
calculations associated with the particle track- 
ing, the interpolation of the gas-phase proper- 
ties, and the source-term evaluation were com- 
puted on the processor of the computational grid 
in which a particle is located. 

This strategy leads to an uniform loading during 
integration but leads to excessive message pass- 
ing. 

• Strategy II: 

The Lagrangian particles were assigned to the 
processor of the computational grid where the 
particle is located. Fig. 9b illustrates a sim- 
ple example of the domain decomposition strat- 
egy adopted for the liquid-phase computations 
where the corresponding gas-flow computational 
domain is divided into equal parts between the 
four available PEs. In this strategy, the La- 
grangian particles are assigned to the processor 


of the computational grid where a particle is lo- 
cated. 

This strategy lead to a non-uniform loading dur- 
ing integration but leads to less message passing. 

Our experience has shown that Strategy II 
seems to work well on different computer platforms: 
both massively parallel computers as well as hetero- 
geneous cluster of workstations. So in the present ver- 
sion of the code, we have opted to implement Strategy 
II over Strategy I. 

16 DETAILS OF PARALLEL 
PERFORMANCE 

The details of the combined parallel perfor- 
mance of the CFD, EUPDF, and LSPRAY codes in- 
volving several different cases can be found in Refs. 
[10-15]. Here, we only summarize briefly the parallel- 
performance results for two different cases. One is 
a 3D test case and more details on this case can be 
found in the reference [13]. For this case, the calcu- 
lations were performed on a computational grid com- 
prising of 8430 tetrahedral elements and 100 Monte 
Carlo PDF particles per cell. The computations were 
performed on one of the NASA Ames Research Cen- 
ter’s parallel computer platforms called Turing which 
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is a SGI Origin work-station with 24 PEs (Proces- 
sor Elements). Table 11 summarizes the CPU times 
per cycle taken by the EUPDF, LSPRAY, and CFD 
solvers vs the number of PEs. Both the CFD and 
PDF solvers show good parallel performance with an 
increase in the number of processors but for the spray 
solver it shows reasonable parallel performance. 

Next, we would like to summarize the results 
from [13, 38] showing only the results of spray compu- 
tations. The results are summarized in Table 12. The 
computations were performed on Turing at NASA 
ARC (Ames Research Center) - it is a SGI Origin 
work-station with a maximum of 24 PEs (Processor 
Elements). The computations made use of an un- 
structured grid with a mesh size of 3600 tetrahedral 
elements. And it made use of about 2,695 Lagrangian 
particles for the spray computations and one hundred 
Monte Carlo particles per element for the PDF com- 
putations. In a given cycle of one global time step, 
dtgi, the spray equations were advanced over one hun- 
dred time steps as given by dt g i/dt m i = 100. The 
first row of Table 12 summarizes the CPU times per 
PE per cycle taken by the spray solver vs the num- 
ber of PEs. As expected, the CPU time goes down 
with an increase in the number of processors. In an 
ideal case, one would expect an inverse reduction in 
cpu time with an increase in the number of proces- 
sors. Here, we don’t get such an ideal performance 
because of the resulting non-uniform distribution of 
spray particles, between various participating proces- 
sors, from the implementation of Stratergy II. To get 
an idea of the spray particle distribution, we have tab- 
ulated the maximum and minimum number of parti- 
cles found between various processors. When we go 
from 1 to 2 PEs, 2097 particles are assigned to one 
and the rest to the second. With four they are dis- 
tributed between 1165 and 118, with eight between 
623 and 14, and with sixteen from 312 to 0. The 
results clearly show that the reduction in the CPU 
time varies almost linearly with the reduction in the 
number of maximum particles. 

17 SUMMARY OF SOME VALIDATION 

CASES INVOLVING BOTH REACTING 
AND NON-REACTING SPRAY 
COMPUTATIONS 

A total of the following nine cases were vali- 
dated: 

1. A reacting methanol spray with no-swirl, 


2. A non-reacting methanol spray with no-swirl, 

3. A confined swirl-stabilized n-heptane reacting 
spray, 

4. An unconfined swirl-stabilized n-heptane react- 
ing spray, 

5. A confined swirl-stabilized kerosene reacting 
spray, 

6. A flashing jet generated by the sudden release of 
pressurized R134A from cylindrical nozzle, 

7. A liquid jet atomizing in a subsonic cross flow, 

8. A Parker-Hannifin pressure swirl atomizer, & 

9. A single-element LDI (Lean Direct Injector) 
combustor experiment. 

The experimental data for the first two cases 
was provided by McDonell & Samuelsen from the 
University of California at Irvine [38] . Both the cases 
are without swirl; one is a reacting case and the other 
is non-reacting. The data for the third and fourth 
cases was provided by Bulzan from the NASA Glenn 
Research Center [39-40]. Both the cases are swirl- 
stabilized reacting cases, one is an unconfined flame 
and the other is confined. The data for the fifth case 
was provided by El Banhawy & Whitelaw from Im- 
perial College [4]. It is a confined swirl-stabilized 
kerosene spray flame. Further details of the last four 
cases can be found in [78]. These cases were chosen 
because of their importance in some aerospace appli- 
cations. 

Here, we provide a brief summary of our val- 
idation effort but a detailed presentation of the re- 
sults and discussion can be found elsewhere in the 
papers [10-13, 78]. The validation is based on some 
3D and axisymmetric calculations involving both re- 
acting and non-reacting sprays. Some of them were 
performed on unstructured grids and the others on 
structured grids. 

The comparisons involved both gas and drop ve- 
locities, drop size distributions, drop spreading rates, 
and gas temperatures. In general, the predicted re- 
sults provide reasonable agreement for both mean 
droplet sizes (U 32 ) and average droplet velocities but 
mostly underestimate the droplets sizes in the inner 
radial region of a cylindrical jet. 

The comparisons also involved the results ob- 
tained from the use of the Monte Carlo PDF method 
as well as those obtained from a conventional CFD 
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solution without the Monte Carlo PDF method. For 
the first case of McDonell & Samuelsen’s reacting 
spray flame, the detailed comparisons clearly high- 
lighted the importance of chemistry /turbulence in- 
teractions in the modeling of reacting sprays [13]. 
The results from the PDF and non-PDF methods 
were found to be markedly different with the PDF 
solution providing a better approximation to the re- 
ported experimental data. The PDF solution showed 
that most of the combustion occurred in a predom- 
inantly diffusion-type of flame environment and the 
rest occurring in a predominantly premixed-type of 
flame environment. However, the non-PDF predic- 
tions showed incorrectly that most of the combustion 
occurred in a predominantly vaporization-controlled 
regime. The Monte Carlo temperature distribution 
showed that the functional form of the PDF for the 
temperature fluctuations varied substantially from 
point to point. The results brought to the fore some 
of the deficiencies associated with the use of assumed- 
shape PDF methods in spray computations. 

18 CONCLUDING REMARKS 


ability to capture the overall structure of a spray 
flame. 

• The source code of LSPRAY-IV will be available 
with NCC as a complete package. 


• This manual provides a complete description of 
LSPRAY-IV - a Lagrangian spray solver devel- 
oped for application with parallel computing and 
unstructured grids. 

• It facilitates the calculation of the multi- 
component liquid sprays with variable properties 
valid over a wide range of low pressure condi- 
tions. 

• It provides the user with a basic understanding 
of the the spray formulation and the LSPRAY- 
IV code structure, and complete details on how 
to couple the spray code to any other flow code. 

• The basic structure adopted for the grid repre- 
sentation and parallelization for the gas side of 
the flow computations follows the guidelines es- 
tablished for NCC. 

• Also, we have extended the joint scalar Monte 
Carlo PDF method to two-phase flows and, 
thereby, demonstrating the importance of chem- 
istry/turbulence interactions in the modeling of 
reacting sprays. 

• Based on the validation studies involving several 
confined and unconfined spray flames, the results 
were found to be encouraging in terms of their 
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APPENDIX 

DETAILS OF THE PRIMARY 
ATOMIZATION MODELS 

1. Sheet Breakup Primary Atomization 
Model (Likely Application: Pressure Swirl 
Atomizer) 

Here, we summarize the details of the sheet 
breakup model taken from Schmidt et al [48]. The 
growth of an infinitesimal disturbance as given by 

p = r; 0 exp(ikx + cot) (1) 

was analyzed based on a linear stability analysis of a 
two-dimensional, viscous, incompressible liquid sheet 
of thickness 2 h which moves through an inviscid, in- 
compressible gas medium. This was analyzed in a co- 
ordinate system moving with the sheet with a relative 
velocity of U where p 0 is the initial wave amplitude, 
k (= 27t/A) is the wave number, and to = to r + icot 
is the complex growth rate. The most unstable dis- 
turbance responsible for the sheet breakup is denoted 
by ft. 

Based on the linearized liquid and gas continu- 
ity and momentum equations subject to the linearized 
boundary conditions at the gas and liquid interfaces, 
a sinuous mode dispersion relation was obtained by 
[49], 


cj 2 [tanh(kh) + Q] + cj [ 4 i 'ik 2 tanh(kh) + 2 iQkU]-\- 
4 vf k 4 tanh(kh) — 4 v 2 k 3 l tanh(lk ) — QU 2 k 2 + = 0 ( 2 ) 

where Q = p g /pi , l 2 = k 2 + uj/vi, and U is the 
relative velocity between the liquid and gas. Inviscid 
analysis also indicates that for low gas Weber num- 
ber (We = p g hU 2 /a) flows, long waves tend to grow 
leading to liquid sheet breakup but for higher Weber 
numbers, short waves produce a maximum growth 
rate followed by breakup. The critical Weber number 
that leads to the transition from the long wavelength 
regime to the short wavelength regime was shown to 
be We g = 27/16. For most modern fuel injection 
systems, the film Weber number is well above this 
critical limit. The growth rate for the sinuous mode, 
ay, based on an order of magnitude analysis of the 
dispersion relation yields, 

2uik 2 tanh{kK) 

ClJ'p 1 — 

tanh(kh) + Q 


4 vf 


k 4 tanh 2 (kh) — Q 2 U 2 k 2 — [tanh(kh) + Q]{—QU 2 k 2 + <rk 3 / pi) 
tanh(kh) + Q 


(3) 


For short waves in the limit of Q « 1 for high-speed 
sheets, it yields 


ay 


-2utk 2 + \Uvfk A + QU 2 k 2 


ak 3 

Pi 


( 4 ) 


Following Dombrowski & Johns [51], the sheet 
disintegration leads to the formation of ligaments 
once the unstable waves reach a critical amplitude 
and Eq. (4) shows that the growth rate of short 
waves is independent of the sheet thickness. The cor- 
responding breakup time r and the breakup length L 
are given by: 


p b = p 0 exp(flr ) => t = \ -In — (5) 

12 p 0 

L = Vt = TjMj) (6) 

i 6 I Iq 

where ft is the maximum growth rate as determined 
by Eq. (4), the term ln( — ) has an assigned value of 
12 as suggested by Dombrowski & Hooper [50], and 
V is the absolute velocity of the liquid. 

The initial diameter of the ligaments is derived 
from a mass balance relationship. For long waves, it 
is assumed that the ligaments are formed from tears 
in the sheet once per wavelength and the resulting 
diameter is given by, 


d-L = 



( 7 ) 


where K sb is the wave number corresponding to the 
maximum growth rate ft as obtained from Eq. (3) 
and the film thickness, h , is calculated from the 
breakup length, L , the radial distance from the cen- 
terline to the mid-line of the liquid sheet at the at- 
omizer exit, r 0 , and the spray angle, 0, as follows: 
h = 2 T p,v(r 0 +L«n( 8 / 2 )) - For short waves, the lig- 
ament diameter is independent of the liquid sheet 
thickness and is assumed to be proportional to the 
wave length associated with the maximum growth 
rate 0 as follows: cIl = 2 j/ c / - where the ligament 
constant, Cl, is equal to 0.5. 

For both long and short waves, Dombrowski & 
Johns [51] developed a linear stability analysis for 
the further breakup from ligaments to droplets based 
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on the Weber’s analysis on capillary instability. The 
analysis shows that the breakup occurs when the 
amplitude of the unstable waves nears the radius of 
the ligament. And the corresponding most unstable 
wavenumber, A'l&, is given by: 


Ku,d,L 


1 + 1-1 

2 2 pi<rd L ) 


(8) 


This analysis thus yields the most probable droplet 
size based on a simple mass balance calculation where 
dl = 3Trd 2 L /K Lb . 


1.1 Application to pressure swirl atomization 


For the pressure swirl atomizer, the initial in- 
jector exit velocity and liquid sheet thickness are cal- 
culated following the approach of Schmidt et al [48] . 
It assumes a uniform velocity profile for the initial 
liquid velocity, V, as given by, 


V = 


max{ 0.7, 


4m 

7T d 2 picos6 



(9) 


where m and 6 are the measured mass flow rate and 
spray angle, respectively, d 0 is the injector exit di- 
ameter, and A p is the pressure drop in the injector. 
Once V is known, the corresponding axial component 
of the sheet velocity is calculated via u = V cosd. 
And the initial sheet thickness h Q is calculated from 
the conservation of mass: 


m = npiuh 0 (d 0 - h 0 ) (10) 

At the point of primary breakup, the actual drop size 
is chosen from a Rosin- Rammler distribution with the 
mean size as given by d b of the sheet breakup model. 
Further movement of the droplets is tracked by mak- 
ing use of a Lagrangian formulation. 

2. Blob Jet Primary Atomization Model 
(Likely Application: Single-Orifice Nozzles) 

Here we summarize the details of the blob 
jet primary atomization model taken from Reitz & 
Bracco [41] and Reitz [52,45]. It applies for a cylin- 
drical liquid jet issuing into an incompressible gas. 
The following dispersion relation was obtained based 
on the stability analysis of a cylindrical liquid surface 
subjected to linear perturbations: 


2 o ,2 r J i( fca ) 2kl h{ka) I[(la) 1 

W 1 ^ I 0 {ka) fc 2 + l 2 I 0 {ka ) I 0 (la ) 1 


ak 

Pia 2 


(1 -fcV){ 


l 2 — k 2 I\{ka ) 
l 2 + k 2 h 0 (ka) 




( 11 ) 


where Jo, I\, and Kq, K\ are the modified Bessel 
functions of the first and the second kinds. 

Reitz [52,45] generated curve-fits of numerical 
solutions to Eq. (11) for the maximum growth rate 
(cj = fl) and the corresponding wavelength (A = A): 


(1 + 0.45Z°' 5 )(1 + 0.4T 0 ' 7 ) , 

9 '° (1 + 0.87Wei-67)0.6 (12) 

0.34 + 0.38 Wef 5 

(13) 

(1 + Z)(l + 1.4T 0 - 6 ) ’ 

where Z = T = ZWe° g 5 , Wei = We g = 

pg ^ a , and Re i = A core region is predicted with 
the blob injection method because there is a region of 
large discrete liquid parcels near the nozzle. Based on 
the jet stability theory, new drops are formed from a 
parent drop or blob. It is assumed that small droplets 
(with radius, r) are formed from the parent drops 
(with radius, a) with drop size proportional to the 
wavelength of the fastest-growing or most-unstable 
wave, 

r = B 0 A {if B 0 A < a) or 
r = min{{3Tra 2 U /2fl) 0 ' 33 , (3a 2 A/4) 0 ' 33 } 

{if B 0 A > a, one time only) (14) 

where B a = 0.61 according to Reitz [52,45]. In the 
above, it is assumed for the {B a A < a) condition 
that small droplets are formed with the dropsize pro- 
portional to the wavelength of the fastest growing 
mode and the second {B 0 A > a) condition applies to 
drops larger than the jet and it assumes that the jet 
disturbance has frequency Q/2i r ( a drop is formed 
each wave period) or that the drop size is determined 
from the volume of liquid contained under one sur- 
face wave. And the rate of change of droplet radius 
due to breakup is given by 


n {^-} 0 - 5 


da/dt = —{a — r)/r (r < a) (15) 

where r is the breakup time, r = 3.726-Bia/Afl, and 
the value for the breakup time constant, B i, depends 
on the injector characteristics and its value ranges 
between 1.732 to 40. And aft = t 0 ) = a 0 is the 
initial drop radius at time, t Q . 
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After the breakup, a new parcel containing 
product drops of size, r, is created and added to the 
computations. This was done as long as the mass of 
the liquid removed from the parent (/9/47r(a 3 — a 3 )/3) 
reached or exceeded 3% of the average injected parcel 
mass and if the number of product drops is greater 
than or equal to the number of parent drops [45]. 
While waiting for sufficient product drops to accu- 
mulate, the parent drop number was adjusted so that 
Na 3 = _/V 0 a 3 but the parent drop number, N a , was 
then restored following the creation of the new prod- 
uct parcel [45]. 

In the case of (B a A > a), the parent parcel 
was replaced by a new parcel containing drops with 
size given by Eq. 14 after a time equal to r (with 
N = N Q a 3 /r 3 ) [45]. This breakup procedure was 
allowed only once for each injected parcel [45]. 

Validation of the model for a single hole orifice 
in a typical diesel engine was demonstrated by Reitz 
and Diwaker [46] and Reitz [52,45]. 


Here we summarize the details of the BLS at- 
omization model taken from Khosla and Crocker [47] . 
In this model both surface shear breakup and col- 
umn breakup modes are included. Before column 
breakup, fragments may be formed due to boundary 
layer stripping depending on the local Weber num- 
ber and q (= piuf / p g v? g ). When the jet reaches the 
column breakup time, the entire jet breaks into frag- 
ments. It also allows for further breakup of the frag- 
ments based on a modified boundary layer stripping. 
And it is followed by a final breakup step based on the 
Rayleigh-Taylor secondary droplet breakup model. 

Fragments are stripped from the liquid column 
if the We g satisfies the following criteria: 

We g > 50Re g /2 9 " 1/0 - 81 and 

We g > 15 (18) 

where 


3. Air Blast Primary Atomization Model 
(Likely Application: Air Blast Atomizers) 

The air blast atomization model is essentially 
based on the idea of pressure-swirl atomization model 
(Section 1.1) as the primary atomization of an air 
blast injector is based on the aerodynamic analysis 
involving the Kelvin-Helmholtz instability of a liq- 
uid jet in an incompressible gas. It, however, differs 
from the pressure-swirl atomization model in the de- 
termination of the initial sheet velocity and thickness, 
which are given by 


V sheet = aVi + (l-a)V g (16) 


S 


r [l 



rh 


7T r 2 pV sheet 


(17) 


where a has a value of 0.12 to 1 depending on the 
fuel filmer characteristics, r is the radius of the fuel 
filmer, and rh is the fuel flow rate. The continuous 
annular sheet is represented by a finite number of 
point injectors located randomly along the circular 
ring of the liquid sheet. 

4. Modified BLS Primary Atomization 
Model (Likely Application: Liquid Jet in a 
Cross Flow) 


Re g — PgdjUg l p g 

we g = u 2 g djp g /ai 

where Re g and We g are based on the gas velocity 
component normal to the liquid jet direction instead 
of the relative velocity. If column stripping does oc- 
cur, the amount of mass removed from the column is 
given by 


where, 


h 


7 rd 


M shed = -Trdpi—UreiAa\/ — At 


A = 


a = 


r Pjl~\ V 3 r Asi V 3 
PZ Mz 

8pi 1 1/2 


r I 

L 3 Au re lpR 


t* = 


Hf. 


(19) 


(20) 

(21) 


(22) 


where is the liquid column drop lifetime. The ad- 
dition of the factor t*/tb causes the shedding rate to 
increase essentially linearly with distance away from 
the injection location. This accounts for the lack of 
shedding close to the injection location and the sub- 
sequent buildup of shedding over the life of the liquid 
column. The shed drop SMD (Sauter Mean Diame- 
ter) is given by 
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SMD = 3A-d 1/2 [— ] 1/4 [_^_] 1/2 (23) 

t* Pa Ur el Pi 

The t b /t* factor is included with the effect of 
producing smaller drops near the injection location. 
However, t b /t* is limited to a minimum value of 2.5 
which is never exceeded for some cases. The amount 
of mass shed is tracked and 10 new parcels are cre- 
ated when the cumulative shed mass after a time step 
exceeds 1% of the mass of the parcel. Each parcel is 
allocated an equal amount of the shed mass and the 
size for each new parcel is selected randomly from a 
uniform distribution between 0.4 and 1.6 times the 
mass mean diameter, MMD. And the drop velocities 
were given by 

u d = u P + 0.3(RND — 0.2) (u g — u p ) (24) 


Vd = v p + +0.3(RND — 0.2)(v g - v p ) (25) 


w d = w p + 0.25 (RND — 0.5) (u re i — w p ) (26) 


the same as Eqs. (24) and (25). The lateral velocity 
component is given by 

Wd = w p + 0A(RND — 0.5 )(u re i — w p ) (28) 

Fragments are further broken into small drops 
according to a modified version of the boundary layer 
stripping model based on the following criteria, 

We g > V Re and 

We g > 15 (29) 

where 

we g = u 2 rel d d p g /<Ji 

Re g = PgddU re i/ p g 

The criteria are generally the same as for column 
stripping except that the dependence on q is not 
needed. Also note that We and Re are now deter- 
mined using the relative velocity instead of the cross 
flow velocity. Again, the mass shed from a fragment 
in a time step and the SMD are given by 


where RND is a random number. Column stripping 
occurs, assuming the above criteria are met, until the 
column breakup time is exceeded. Parcels created 
through the column stripping mechanism are con- 
sidered to be fragments which may undergo further 
breakup as discussed below. First, though, the col- 
umn breakup mechanism, which also produces frag- 
ments, is described. The liquid jet column breakup 
time is given by 

t b = A b We 0 62 t* (27) 

The constant, A b , has a value of 25. Since the present 
model includes a fragment breakup step after the col- 
umn breakup, the We dependence for the breakup 
onset was retained. 

After the column breakup time is reached, the 
column is broken into 18 new parcels with MMD = 
0.45dj. The new parcels are also designated as frag- 
ments. The size of the fragments still tend to be 
large, so the ultimate drop size from the primary 
breakup process is mostly determined from the frag- 
ment breakup process. The size distribution is the 
same as described above for the column stripping. 
The cross flow and normal velocity components are 


M shed = 1.2ndpiu re iAay —A t (30) 

SMD = 3.6d 1/2 [—] 1/4 [— ^ — 1 1/2 (31) 

where A and a are given by Eqs. (20) and (21), re- 
spectively. The new droplet velocities are given by 
Eqs. (24), (25), and (28). The broken fragments pro- 
duce 3 new parcels with size distribution the same as 
described above when the shed mass from the frag- 
ment exceeds 20% of the fragment mass. A fragment 
can continue to breakup until it no longer meets the 
criteria of Eq. (29) or until its size is lower than the 
newly created drops. Once the fragment breakup pro- 
cess is complete, drops may breakup further based on 
a Rayleigh-Taylor secondary breakup method. 

Khosla and Crocker [47] applied the model to 
predict the properties of Jet A-l kerosene fuel in- 
jected into a cross-flowing air stream. 

DETAILS OF THE SECONDARY 
DROPLET BREAKUP MODELS 
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Recent studies by Reitz et al [55, 52] have ex- 
amined the breakup of single droplets moving in 
a transverse, high-velocity air jet. The high-speed 
photography provided new insights into the details 
of the breakup mechanism of a single drop. The 
droplet breakup regimes are classified as bag, strip- 
ping (shear) and catastrophic (surface wave) based 
on an increasing size of Weber number. In the bag 
breakup mode (at low Weber number), the drop is 
flattened by the aerodynamic pressure, then turned 
inside out, forming a thin hollow bag which is tied to- 
gether with a circular belt-like structure on the wind- 
ward side. The bag eventually bursts into smaller 
liquid fragments, whereas the belt decays into larger 
ligaments and droplets. In the stripping regime thin 
sheets or ligaments of fluid are continuously shed from 
the periphery of the distorting parent drop as a con- 
sequence of a K-H instability, causing these sheets to 
disintegrate into tiny droplets. This process always 
leaves a coherent residual parent drop. The catas- 
trophic drop breakup takes place in two stages lead- 
ing to a collection of larger and tiny product droplets: 
Large amplitude long-wavelength waves caused by 
drop deceleration induce a R-T instability on the 
flattened drop which leads to a breakup into large 
product droplets, while at the same time short sur- 
face waves induce a K-H instability on the windward 
side of the parent drop resulting in a collection of 
much smaller product droplets. In diesel or other 
high-pressure gas-turbine sprays the droplets span a 
wide range of velocities and hence Weber numbers, 
and thus it is expected that all three droplet breakup 
mechanisms are relevant in the breakup modeling. 

In what follows, we provide some details of 
the secondary droplet breakup models contained in 
the CFDRC/UW atomization module: (1) Rayleigh- 
Taylor, (2) TAB, and (3) ETAB. These details are 
taken from [52-56, 42-44]. 

1. Rayleigh- Taylor Secondary Droplet 
Breakup Model 

Here we summarize the Rayleigh-Taylor sec- 
ondary breakup model developed by Patterson and 
Reitz [57]. It is based on the analysis developed 
by Taylor [53,57] that accounts for the disturbances 
caused by droplet deceleration. In the Rayleigh- 
Taylor breakup mechanism, the breakup wavelength, 
A, is given by 

A = Q.2-K^Za/\ud\{pi - p g ) (32) 


where ii,j is the drop deceleration (= | Cd ^ a J 7 , Cd is 
the aerodynamic drag coefficient, U is the drop rela- 
tive velocity, & d is the drop diameter) . Furthermore, 
the breakup wavelength A is limited by 


A = max(0.8d, A)[l + 0.2(RND — 0.5)] (33) 

and the breakup time, tb,RTi is given by 


tb,RT = Cfreq > 0.5^(pi + Pg)( 


\Ud\(pi~Pg) 


1.5 


(34) 

However, the value assigned for the constant, Cf req , 
depends on the droplet classification - parent, prod- 
uct, or default. For more details, one can refer to 
Patterson and Reitz [21]. After the breakup, no new 
drop parcels are created and there is no change in ve- 
locities between the parent and product drops. How- 
ever, the drop number in a given parcel changes to 
n product as given by n paren t{d parent /A) a due to the 
change in the sizes between the parent and product 
drops. 


2. The TAB Secondary Droplet Breakup 
Model (Likely Application: 

Lenticular-Shaped Droplet Deformations) 

In an attempt to provide a description of the 
droplet and jet disintegration in the modeling of 
diesel sprays, O’Rourke & Amsden introduced the 
TAB model [42] . Here we summarize the TAB model 
taken from [43-44]. 

The TAB breakup model is based on Taylor’s 
analogy between an oscillating, distorting drop and a 
spring-mass system. A detailed analysis of this model 
together with a discussion of its numerical implemen- 
tation can be found in [42 & 58]. In this model, the 
drop motion is governed by a linear ordinary differ- 
ential equation for a forced, damped harmonic oscil- 
lator. The forcing term is given by the aerodynamic 
droplet-gas interaction, the damping is due to the 
liquid viscosity and the restoring force is supplied by 
the surface tension. The parameters and constants 
have been determined partly from theoretical consid- 
erations and partly from experimental observations. 

The TAB model describes the distortion of the 
drop by the deformation parameter, y = 2x/a, 

where x denotes the increase in the radius increase 
of the equator from its equilibrium position and a 
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is the drop radius. The equation for the distortion 
parameter y is given by 


V + 


5 m 

Pia 2 


y + 



2p g \u\ 2 

3 pia 2 


(35) 


Assuming a constant relative drop-gas velocity, 
U (which is satisfied in the numerical solution process 
during a given time step), the solution to Eq. (35) is 
given by 


. , We 

M = 12 + 


=i n We-, 

e *d ( [y( 0 ) — J cosojt 


-y 

L u; 


2/(0) - if - 

utd 


smio 


l ) 


(36) 


Here we provide some details of the ETAB 
model taken from Tanner [43-44]. The ETAB model 
is based on the following modifications to the stan- 
dard TAB model: (1) the droplet disintegration is 
modeled via an exponential law which relates the 
mean product droplet size to the breakup time of the 
parent drop; and (2) an energy balance consideration 
between the parent and product droplets yields an 
expression for the normal velocity component of the 
product droplet. 

When the breakup condition of We > 
W e cr it = 6 and y(t) = 1 is met, then the parent 
drop breaks up into a collection of product droplets, 
subject to a size distribution function which, in gen- 
eral, depends on the breakup mechanism. In the 
ETAB model, the rate of product droplet generation, 
dn(t)/dt, is given by 


where Ife = p g aU 2 /a and 


td 


2pia 2 
5 pi 


(37) 


8a 1 

Pia 3 t 2 


(38) 


In this model, it is assumed that a necessary 
condition for drop breakup is reached when We > 
We cr it . And the value for the critical Weber num- 
ber is determined experimentally to be 6. For an 
inviscid liquid with initial deformation conditions 
2/(0) = 2/(0) = 0, the solution to Eq. (35) 

reduces to y(t) = We( 1 — cosw 0 t)/ 12, where 

co 2 = [44]. Consequently, the breakup occurs 

when y(t) > 1. The drop size after breakup is de- 
termined by an energy balance equation between the 
parent and the product droplets which equates the 
surface energies of parent drops with the energies of 
product drops due to oscillation and distortion. Also, 
the product droplets are initially equipped with the 
additional deformation velocity x = ay/2, which 

acts normal to the path of the parent drop and is 
responsible for the formation of the spray angle. 

One major advantage of this model is that it is 
based on a simple linear equation and it can be used 
effectively to describe the lenticular-shaped droplet 
deformations as observed in the experiments of [56 & 
54], 


3. The ETAB Secondary Droplet Breakup 
Model (Likely Application: High-Pressure 
Diesel Engine) 




(39) 


where n{t) = mo/rh(t ) and mo is the mass of the 
parent drop and to the mean mass of the product 
droplet distribution. Utilizing the fact that dn/dt = 
— ( mo /in 2 ) (dm /dt ), leads to the breakup law which 
relates the product drop size to the breakup time as 
determined by the TAB model. 


— = -3K b m (40) 

The breakup constant K b depends on the 
breakup regime and is given by parent drop prop- 
erties only. Bag breakup occurs if We = Wet and 
stripping breakup happens if We > We t . And it is 
given by 


I\ b = k\LO if We < We t or 

K b = k 2 uVWe if We > We t (41) 

The values for Wet, fci and k 2 have been determined 
experimentally and has been set to k\ « fc 2 = 1/4.5 
and We t = 80. 

In this model, a uniform product droplet size 
distribution is assumed. It is also noted that the 
choice of uniform distribution is not expected to be 
realistic but may produce good approximations when 
averaged over many drop breakups, because parent 
drops of different sizes and Weber numbers will in 
general yield a wide range of duct droplet sizes. With 
this assumption, Eq. (40) becomes 

- = e~ Kbt (42) 

a 


NASA/CR— 2012-217294 


50 



where a and r are the radii of the parent and product 
drops, respectively. 

After breakup of a parent drop, the initial de- 
formation parameters of the product droplets are set 
to y( 0) = y(0) = 0. Also, the product droplets are 
initially supplied with a velocity component perpen- 
dicular to the path of the parent drop with a value 
vt = Ax, where A is a constant determined from the 
following energy balance consideration. The energy 
of the parent drop is the sum of the surface tension en- 
ergy and the droplet deformation energy. The second 
is computed as the product of the aerodynamic drag 
and the drop deformation at the stagnation point, 
estimated to be 5a/9. This leads to 


Eparent = 47rcra 2 + 5irc d p g a 3 \U\ 2 /18 (43) 

And the energy of the product droplets in the frame 
of reference of the parent drop is given by 


Eproduct = 47 Taa 3 /r S MR + A 2 npia 5 y 2 /6 (44) 

where the Sauter mean radius, tsmr , enters via the 
relation r 2 = a 3 /tsmr- From Eqs. (43) and (44) 
one obtains the relation 


A 2 = 3[1 - a/rsMR + 5c d We/72}u 2 /y 2 (45) 

where to 2 = 

Pia s 

Tanner [43-44] analysis yields an approximate 
value of 0.69 for A showing that only 70% of the par- 
ent drop deformation velocity goes into the normal 
velocity component of the product droplets, where as 
it is 100% in the standard TAB model. 

Also, the characteristic time, r (= A~), for 
breakup in Eq. (42) for an inviscicl liquid (pi = 

0) is given by a± yj if We < We t or 
if We > We t , where the suggested val- 
ues are for a± = (VSki)^ 1 and a.i = (v / 8/c2)^ 1 - 

The application of the ETAB model in the sim- 
ulation of a high pressure liquid jet breakup can be 
found in [43-44]. 
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